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Non-Propagating Cracks in Vee-Notched 
Specimens Subject to Fatigue Loading 


N. E. FROST 


(Mechanical Engineering Research Laboratory, East Kilbride, Glasgow) 


Summary: Reversed direct stress and rotating bending fatigue tests have been 
carried out on Vee-notched specimens of aluminium alloy, nickel chromium 
steel and mild steel. 


Diagrams are presented showing the relationship between the geometric 
stress concentration factor K, and the strength reduction factor K;. It was found 
that non-propagating cracks formed in the roots of the sharper notches. These 
cracks formed at or above some critical value of K,, the value depending on the 
specimen material. Below the critical value of K,, cracks did not form unless 
the applied nominal stress exceeded that at the fatigue limit, and a crack, once 
formed, continued to propagate until the specimen failed. Above the critical 
value of K,, non-propagating cracks formed at nominal stresses less than the 
fatigue limit, the stress having to be increased to the latter value in order to 
propagate the crack. This critical value of K, coincided with the maximum K, 
value realised. It would appear that K, equals K, up to a certain value of 
K,; there is then a transition where K, reaches a maximum at the critical value 
of K,. Increasing the value of K, above the critical value causes no further 
increase and may tend to decrease the value of K;. 


The conclusions drawn apply only when the fatigue stresses are completely 
reversed, i.e. the mean load is Zero. 


1. Introduction 


The problem of assessing the effect of a stress raiser on the fatigue properties 
of a material has been the subject of a considerable amount of work. Simply stated, 
the problem is to derive the value of K; (the strength reduction factor represented by 
actual behaviour) from the value of K, (the geometric stress concentration factor) 
for the stress raiser in question. K, is defined as the ratio of the maximum stress in 
the neighbourhood of the stress raiser to the nominal stress assuming perfect 
elasticity and can be obtained in certain cases by mathematical analysis or can be 
determined experimentally, usually by photoelastic techniques. Under conditions of 
fatigue loading, maximum stress values less than the yield point of the material can 
cause failure, and hence the value of K, should be important. In static loading, on 
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Fatigue specimens for direct stress machines (with notch). 


the other hand, plastic flow will occur at stresses above the yield point and will alter 
the elastic stress field around the stress raiser before final failure occurs. As the value 
of K, is based on the maximum principal stress, reference to this factor assumes that 
failure depends on the maximum principal stress. However, combined stress 
fatigue tests show that for ductile materials failure obeys the von Mises criterion 
more closely than that of maximum principal stress. Peterson‘, assuming the 
von Mises criterion of failure to apply, has computed values of a geometric elastic 
stress concentration factor, which he terms K,’, for several types of stress raisers. 
The value of K,’ is never greater than K,, and the maximum possible difference 
between the two factors is about 15 per cent. As differences much larger than this 
are found between K, and K;, K, is usually preferred to K,’ on the grounds of 
greater simplicity. 


Much work has been devoted to the experimental determination of K;, but 
mostly in small rotating beam machines, and comparatively few tests have been 
made in reversed direct stress using large specimens. 


It is not intended to go into detail regarding the relationships that have been 
suggested to connect K, and K; (e.g. the conception of a notch sensitivity index 
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Fig. 2. 
Fatigue specimen for rotating bending fatigue machine (with notch). 


q=(K;-— 1)/(K,— 1) or Kuhn and Hardrath’s work” based on Neuber’s™ elementary 
structural unit); examples of these have been collected and discussed by Dolan and 
Yen’. Several Authors, e.g. Phillips and Heywood, Hikata® and Dolan and 
Hanley™, have found that, for geometrically similar test pieces, K; increases with 
absolute size of test piece. This is termed “size effect.” In contrast to these 
findings, Hyler, Lewis and Grover“ and Ros and Eichinger®’, appear to find no 
appreciable size effect for the normal range of sizes met in practice. However, Ros 
and Eichinger™ state that a size effect is apparent for very sharp discontinuities. 


A further point which will be brought out later is that several authors, 
e.g. Mann"), Hyler, Lewis and Grover“, and Ros and Eichinger®, have found that 
the sharpest notch does not necessarily have the greatest weakening effect. Mann’, 
working with an aluminium alloy 24 S-T in rotating beam tests, found that a notch 
with a K, value of 5-8 gave a maximum K, value of 3:2, K; decreasing with further 
increase in K,. Frost and Phillips'''’, from work on the fatigue strength of cracked 
specimens, suggested that the sharpest notch is not necessarily the most damaging. 
In the present paper the effect on the relationship between K, and K;, of the forma- 
tion of a non-propagating crack at the root of the stress raiser is considered. The 
presence of such non-propagating cracks has been reported by Fenner, Owen and 
Phillips’” for mild steel and Frost''*’ for aluminium alloy. For both materials 
cracks were found at the roots of Vee-notched specimens at stresses from about 
‘++ 1 ton/in.? up to the nominal fatigue limit in specimens run for 30 - 100 x 10° cycles 
and which normally would be considered unbroken. Further experimental data 
are presented in the paper regarding the relation of the notch dimensions to the 
formation of such non-propagating cracks. 


NOTATION 
K, strength reduction factor 
K. geometric elastic stress concentration factor 
K,’ geometric elastic stress concentration factor (Ref. 1) 
a_ radius of minimum cross section of specimen 
notch depth See Fig 3 
r root radius 
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2. Experimental Details 


Three materials were used for the investigation: —a 0-15 per cent carbon mild 
steel, a nickel chromium steel (EN 26) and an aluminium alloy (BS L65). Details 
of their chemical composition and mechanical properties are given in Table I. 


Fatigue tests on notched specimens were carried out under both reversed direct 
stress and rotating bending. For both these types of loading the mean load of the 
cycle is zero and hence the applied nominal stress cycle is completely reversed. The 
types of specimen used in direct stress are shown in Fig. 1, while that used in 
rotating bending is shown in Fig. 2. Fig. 3 illustrates the dimensions a, ¢ and r used 
to describe the notch. In all cases a Vee type of notch with an included angle of 
55° was used. Great care was taken in the machining of the notches; the accuracy 
of the root radii obtained can be seen from the various sectional photographs given 
in the paper. 


The reversed direct stress tests were done in a Schenck type fatigue testing 
machine running at approximately 1,800 cycles/min., while the rotating bending 
tests were made in a 10,000 in. Ib. capacity Sonntag four-point loading machine 
running at 3,000 r.p.m. 


Specimens which were unbroken after running for 10-30 10° cycles were 
examined microscopically as detailed in Ref. 13. Briefly, the section of the specimen 
containing the notch was sawn across two mutually perpendicular diameters, thus 
exposing four different faces for viewing. The specimens were polished and etched 
in the usual manner. 
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Fig. 3. 


Dimensions of notch in terms of a, t and r. 
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TABLE I 
DETAILS OF MATERIALS USED 


Chemical | H na] Tensile 4 Ol per cent | Elongation 
Material Symbol analysis | treatment | strength | proof stress |on4 (area) 
) | (percent) | | (tons/in.2) | (tons/in.?) (per cent) 
015 C | 
. } 045 Mn | Normalise | 29 22 30 
Mild steel TON 
0-014 920°C (yield stress) 
0-015 P 
0:43 C 
2°64 Ni | 
0:75 Cr 850°C i 
Nickel 0°65 Mn | oil quench : 
chromium RSU 0:58 Mo | 63 54 24 
steel 0-05 Va 650°C | (yield stress) 
0°32 Si air cool 
0:012 P 
0012S 
44Cu 505°C 
Aluminium| MBCG\| | 0:7 Mg water 39 28 8 
alloy MBHQ{ | 0-7Si quench age 
0-6 Mn 160-190°C 
TABLE II 
VEE-NOTCHED ALUMINIUM ALLOY SPECIMENS 
Values of Fatigue Limit, K,, K, 
and Remarks on the Microscopic Examination of Notch Root 
Fatigue 
r a t endurance K K | Microscopic examination and 
(in.) (in.) (in.) limit : , | photograph figure number 
(tons/in.?) 


Reversed Direct Stress 


0-05 0-65 0-2 <3 3-3 | 3-2 No crack detected (Fig. 6). 
0:02 0°65 0-2 2 oy 4:85 | 4:85 | Nocrack detected (Fig. 7). 
0-008 0°65 0-2 7:3 6°5 No crack detected (Fig. 8). 
0-004 0-65 0:2 + 1-5 10 6°5 Small crack inclined at 45° to 

applied load (Fig. 9). 
0-002 0°65 0-2 +25 | 14 3-9 Crack at 90° to applied load 
(Fig. 10). 
0-0005 0°65 0:2 +26 27 3-7 Crack at 90° to applied load. 
(approx.) Multi-cracks formed (Fig. 11). 
Rotating Bending 
0:02 0°55 0-2 No crack detected (Fig. 12). 
(approx.) 

0-008 0:70 0-05 +3°4 49 2:85 No crack detected (Fig. 13). 
0-002 0:55 0-2 Cracks at 45° to applied load 

(Fig. 14). 
0°0002 0°55 0-2 +36 34 Crack at 90° to applied load 
(approx.) (Fig. 15). 
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S-N curve for plain specimens of aluminium alloy. 


Experimental Results* 
3.1. ALUMINIUM ALLoy B.S. L 65 


Preliminary fatigue tests under both rotating bending and reversed direct stress 
were done on plain specimens. The results are shown in Fig. 4. Wo6ohler specimens 


of 0-3 in. diameter were used for the rotating bending tests, and specimens 1-38 in. 
diameter for the reversed direct stress tests. There is little difference between the 
two sets of points and a common curve, giving a fatigue limit for the material of 
+9-7 tons/in.? at an endurance of 50 x 10° cycles, has been drawn. 


Table II gives details of the different types of notches tested under reversed 
direct stress and rotating bending. This table also shows the corresponding estimated 
fatigue limits, the values of the geometric stress concentration factor K, and the 
strength reduction factor K;, and the results of the microscopic examination of the 
material at the root of the notch of the unbroken specimens. A typical S-N curve 
for one particular notch shape is given in Fig. 5. K;, is obtained from Neuber’s 
analysis’, while K; is the ratio of the plain fatigue limit to the notched fatigue limit. 
Photomicrographs of the material at the root of the notch of unbroken specimens for 
each type of notch are given in Figs. 6-15. 


*Full experimental details are given in Ref. 18. 
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Fig. 5. 
S-N curve for Vee-notched aluminium alloy specimens. 
a=0°65 in. t=0°2 in. r=0-008 in. 


Reversed direct stress. 


3.2. MILD STEEL 0:15 PER CENT C 


The results of the preliminary fatigue tests on plain specimens are shown in 
Fig. 16. The conditions were similar to those of the aluminium alloy tests, except 
that the direct stress specimens were 0-977 in. diameter. The difference between 
the results obtained from the two types of loading could be explained on a stress 
re-distribution basis, as considered by Forrest and Tapsell“”. The lower fatigue 
limit of +14-0 tons/in.’, i.e. the fatigue limit for the direct stress tests, has been 
used as the basic fatigue limit for the material. This value agrees quite closely 
with that obtained by Phillips and Heywood"? for similar material and consequently 
their results, where appropriate, have been included with those obtained in 
this paper. 


Table III (p. 13) presents information corresponding to that given in Table II 
for the aluminium alloy. A typical S-N curve for a particular notch shape is shown 
in Fig. 17. 


Photomicrographs of the material at the root of the notch of the unbroken 
specimens are shown in Figs. 18-22. 


3.3. NICKEL CHROMIUM STEEL (EN 26) 


Unfortunately, material for only a few tests on this steel was available. As 
with the mild steel, some results obtained by Phillips and Heywood for a similar 
steel have been included. These authors found little difference between the fatigue 
limits determined by rotating bending and reversed direct stress tests on plain 
specimens, and an average value of +36°5 tons/in.? has been taken from their 
results. 
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Fig. 17. 


S-N curve for Vee-notched mild steel specimens. 


a=0°65 in. 


t=0-2 in. 


r=0-004 in. 


Reversed direct stress. 
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Vee-notched aluminium alloy specimens. 
Reversed direct stress. 


ter 
Fig. 6. Fig. 7. 
Specimen MBCG 2L x 25* Specimen MBCG 2G x 100 
r=0:05 in. r=0-02 in. 
+3 tons/in.? 27°32 x 10° cycles unbroken. +2 tons/in.? 31-52 x 10° cycles unbroken. 
Age. 
Fig. 8. Fig. 9. 
Specimen MBCG IIN x 100 Specimen MECG 18A 500 
r=0-008 in. r=0-004 in. 
+1°5 tons/in.? 33:26 x 10° cycles unbroken. + 1:5 tons/in.? 30 x 10° cycles unbroken. 


*The magnifications quoted in all the microphotographs refer to the originals. !n this printed 
version they are about 20 per cent less. 
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Vee-notched aluminium alloy specimens. 
Reversed direct stress. 


: 

‘ 

¢ 4 
Fig. 10. Fig. 11. 
Specimen MBCG 14M 500 Specimen MBHQ 3C x 500 
r=0-002 in. r=0-0005 in. (approx.) 

+2 tons/in.? 28:05 x 10° cycles unbroken. +2°5 tons/in.? 29°26 x 10° cycles unbroken. 


Rotating bending. 


> 
Fig. 12. Fig. 13. 
Specimen MBCG 7QI1 x 100 Specimen MBHG 16] x 250 
r=0-02 in. r=0-008 in. 
+ 3 tons/in.* 30-3 x 10° cycles unbroken. +3 tons/in.? 17°35 x 10° cycles unbroken. 
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Vee-notched aluminium alloy specimens. 
Rotating bending. 


Fig. 14. Fig. 15. 
Specimen MBCG 7LI x 500 Specimen MBHQ 16F x 500 
r=0-002 in. r=0°0002 in. (approx.) 
+ 2°75 tons/in.? 30-05 x 10® cycles unbroken. + 3-25 tons/in.? 17:16 x 10° cycles unbroken. 


Vee-notched mild steel specimens. 
Reversed direct stress. 


Fig. 18. 
Specimen TON 4F x 200 
r=0-05 in. 
+ 4:25 tons/in.2 10°63 x 10° cycles unbroken. 
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Vee-notched mild steel specimens. 
Reversed direct stress. 


Fig. 19. Fig. 20. 
Specimen TON 4G x 250 Specimen TON 4P x 250 
r=0-01 in. r=0:004 in. 
+ 3°75 tons/in.? 23:39 x 10° cycles unbroken. + 3-25 tons/in.? 20°8 x 10’ cycles unbroken. 


Rotating bending 


Fig. 21. Fig. 22. 
Specimen TON 1Q x 250 r=0-004 in. 
r=0-004 in. Specimen TON 3D ~x 250 
+4 tons/in.? 12°75 x 10’ cycles unbroken. +5 tons/in.? 15°97 x 10’ cycles unbroken. 
12 The Aeronautical Quarterly 


“as. 
+ 


NON-PROPAGATING CRACKS 


TABLE III 
VEE-NOTCHED MILD STEEL SPECIMENS 
Values of Fatigue Limit, K,, K, 
and Remarks on the Microscopic Examination of Notch Root 


r a | t Fatigue | * | Microscopic examination and 
(in.) (in.) |  (in.) limit | photograph figure number 
(tons/in.*) | | 
j Reversed Direct Stress 
0:05 0°65 0:2 |. No cracks (Fig. 18). 
0:025* 0°65 0-2 | 44 | 38 No examination. 
0-01 0°65 0:2 +375 | G6 | 7 Crack at 90° to applied load 
6 | (Fig. 19). 
0:004 0-65 0:2 + 3-4 | 10 41 Crack at 90° to applied load 
(Fig. 20). 
‘ 0-002* 0°65 0:2 144 | 38 Cracks present. See Ref. 12. 
Rotating Bending 
0-004 0°55 0:2 8:0 Crack at 90° to applied load 
(Fig. 21). 
0-004 065 | O1 | +5 TS | 2s Crack at 90° to applied load 
| (Fig. 22 
TABLE IV 


VEE-NOTCHED NICKEL CHROMIUM STEEL SPECIMENS 
Values of Fatigue Limit, K,, K; 
and Remarks on the Microscopic Examination of Notch Root 


r a t Fatigue | Microscopic examination and 
(in.) (in.) (in.) limit » | K, photograph figure number 
(tons /in.?) 
Reversed Direct Stress 
0-005 0°425 0:02 +80 46 4-6 No crack detected. 
0:005 0°425 0:2 +46 8-0 | 19 No crack detected. 
Rotating Bending 

0-004 0°55 0-2 +6°25 so | 58 No crack detected. 

0-010 | 0-55 0:2 +7-25 52 | 50  |No crack detected. 


Table IV presents information for this material corresponding to that given in 
Tables II and III. A typical S—N curve is shown on Fig. 23. 


For the few different notches tested, no non-propagating cracks could be 
detected in any of the unbroken specimens. In view of this, no photomicrographs 
are given for this material. 


4. Discussion of Results 


As most results have been obtained with the aluminium alloy, these will be 
discussed in detail, and the mild steel and nickel chromium steel results dealt with 
more briefly. 


* Results taken from Ref. 5. 
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Fig. 23. 
S-N curve for Vee-notched nickel chromium steel specimens. 
a=0°425 in. t=0°2 in. r=0-005 in. 


Reversed direct stress. 


The main purpose of the experimental work was to study the formation of cracks 
and for this reason a wide variety of notch shapes was used. Although only 
relatively few specimens were tested for each individual notch shape, it is considered 
that the estimated fatigue limits are not seriously in error. 


For the aluminium alloy a series of tests in reversed direct stress was carried 
out with the values of a (radius of minimum cross section) and ¢ (notch depth) kept 
constant and only r (root radius) varying. Fig. 24 shows the fatigue limit plotted 
against root radius for these tests. The point at zero root radius is taken from 
results by Frost and Phillips" relating to cracked specimens. Although this point 
cannot, for several reasons, be taken as strictly comparable with the others, it 
accords with their general trend and certainly represents an upper limit to the true 
value for zero notch radius. It can be seen that the fatigue limit decreases with 
decreasing root radius to a root radius value of about 0-0055 in. and then increases 
again. As mentioned previously in Section 1, the fact that a specimen containing 
a notch of sharp root radius ¢an be stronger than one with a less sharp root radius 
has been noted by several authors’: *:'’’, but no completely satisfactory explanation 
has been put forward. 


Now, considering the vertical line on Fig. 24 through the critical root radius, 
microscopic examination of the material at the root of the notch of unbroken 
specimens (i.e. those run at the fatigue limit or just below) with root radii greater 
than this value (to the right of the line) showed no breakdown of the material, as 
can be seen from Figs. 6, 7 and 8. However, equivalent specimens with smaller 
root radii (to the left of the line) were all found, on examination, to contain non- 
propagating cracks (see Figs. 9, 10, and 11). Hence the increase in fatigue strength 
obtained with root radii to the left of the vertical line is coincident with the 
formation of non-propagating cracks. 


The notch shown in Fig. 11 was intended to be as sharp as possible, i.e. as 
sharp as could be machined. The resulting notch was rather irregular in profile 
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but, as it was obviously sharper than the next larger size, the specimens were tested 
as they were. For the equivalent sharp notch tested in rotating bending, a different 
machining technique was employed, and a much better notch was obtained (Fig. 15). 


The formation of non-propagating cracks in the aluminium alloy material has 
been discussed previously by the author’). For the purpose of this paper it is 
sufficient to state that the cracks formed extended around the whole periphery and 
reached a reasonably uniform depth. With the reversed direct stress tests, this 
uniformity of crack growth is achieved only when the loading is truly axial. 


Further work is required before a detailed explanation can be given of the 
causes and critical conditions necessary for a crack to remain dormant. Any 
explanation must involve consideration of the stress gradient generated by the 
notch and the variation of this stress gradient as the crack forms at the notch root. 
However, as no cracks were detected to the right of the line on Fig. 24, it would 
appear that the nominal stress has to be increased with increasing notch radius in 
order to initiate the crack, which then propagates automatically; while to the left 
of the line a non-propagating crack is formed, and the nominal stress needs to be 
increased with decreasing notch radius in order to cause the crack to propagate and 
the specimen to fail. 


If, in Fig. 24, the point representing the stress necessary to initiate a crack at 
0-002 in. radius"*) is plotted, it lies on a feasible continuation (as shown by the 
dotted line) of the curve to the right of the critical radius back to the origin. The 
difference between this dotted line and the full line to the left of the critical radius 
is the increase in nominal stress required to propagate the crack formed at the stress 
represented by the dotted line. 


VALUE FOR CRACKED SPECIMEN 
FROM REF.) 


FOR ALL NOTCHES a=0-65IN. t*0-2 IN. 


RADI GREATER THAN THIS VALUE 


DO NOT PRODUCE NON-PROPAGATING CRACKS 
> | ALL RADII LESS THAN THIS VALUE 
2 | PRODUCE NON-PROPAGATING CRACKS 
+44 
= | 
5 
w £3} 
5 
MINIMUM NOMINAL FATIGUE LIMIT 


| 


i iL i L i 
0-02 0-05 0-04 0-05 0-06 
ROOT RADIUS (IN) 
Fig. 24. 


Fatigue limit of aluminium alloy Vee-notched specimens in reversed direct stress plotted 
against root radius. For all notches a=0-65 in., =0-2 in. 
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All the experimental data obtained for the aluminium alloy are plotted on 
Fig. 25, which shows the geometric factor K, plotted against the strength reduction 
factor K;. Consider first the direct stress tests only; it can be seen that the points lie 
on the K,=K, line up to a certain value of K,, reach a maximum value of K;, and 
then decrease. All unbroken specimens having a K, value to the right of the 
vertical line drawn through the maximum K; point were found to contain non- 
propagating cracks (Figs. 9, 10 and 11). No cracks could be detected in unbroken 
specimens with the K, values lying to the left of the line (Figs. 6, 7, and 8). Hence 
it would appear from the direct stress tests that the value of K; equals K, up to a 
certain value of K,; then there is a transition where K; reaches a maximum at some 
critical value of K,. Beyond this critical value of K, non-propagating cracks form. 
and the value of K; then decreases with increasing values of K,. 


Considering now the rotating bending tests, it can be seen that although the 
general shape of the curve is similar, the maximum K, value is much smaller than 
for direct stress, although it occurs at about the same critical value of K;. Again, 
unbroken specimens with K, values to the right of the line were found to contain 
non-propagating cracks (see Figs. 14, 15) while no cracks could be detected in 
unbroken specimens with K, values to the left of the line (see Figs. 12, 13). Also 
included in the curve are numerous points which are taken from the literature, the 
numbers adjacent to the points indicating the reference. It appears that all the 
points lie reasonably about the one curve, despite the fact that the tests relate to 
various specimen and notch sizes. It is not unreasonable to assume that the two 
curves tend to become asymptotic to the K; value for a cracked specimen"”’. 
Although it is not suggested that the dotted horizontal line representing the cracked 
specimen result is strictly comparable, it can almost certainly be regarded as a lower 
limiting value. 


It is not meant to imply that the shapes of the two curves on Fig. 25, especially 
to the right of the vertical line through the critical K, value, have been precisely 


NON-PROPAGATING CRACKS 
FORMED IN ALL SPECIMENS 
TO RIGHT OF LINE 


NO CRACKS i 
DETECTED IN | 
SPECIMENS TO | 
LEFTOF.LINE | 
| 
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Fig. 25. 


Relation between K, and K, for aluminium alloy. 
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Fig. 26. 
Relation between K, and K, for mild steel. 


determined. The position of the vertical line in Fig. 25 is reasonably well defined, 
but to the right of this the curves as drawn must be regarded merely as indicative 
of general trends. Investigation of this portion of the curve in detail would involve 
considerable difficulties in specimen preparation, because of the small root 
radii required. 


Figure 26 shows K;, plotted against K, for the mild steel specimens. Again 
two separate curves are obtained, one for reversed direct stress and one for rotating 
bending. The value of K, corresponding to the maximum K, values appears to be 
the same for both methods of test. As for the aluminium alloy, all unbroken 
specimens with a K, value greater than the critical value were found to contain non- 
propagating cracks, whether tested in reversed direct stress (Figs. 19, 20 and 
Ref. 12), or rotating bending (Figs. 21 and 22), and no cracks could be found in 
unbroken specimens with a K;, value less than the critical value (Fig. 18). 


Again the curves to the right of the critical K, value can only be considered as 
showing general trends. For the mild steel it appears that, above the critical value 
of K,, the value of K; may remain constant. It certainly does not decrease as 
markedly as was found with the aluminium alloy. The horizontal line representing 
the K; of the cracked specimen” relates to the direct stress value, and again can 
only be regarded as a lower limiting value. 


Table V gives the critical values of K, obtained for the three materials and the 
corresponding maximum K, values, together with the limiting K;=K, values. 


Although insufficient experimental results are available to plot K, against K; 
for the nickel chromium steel, the results obtained when supplemented by data 
from Refs. 3 and 5 make it possible to give the approximate values included in 
Table V. These values show the differences between the two types of steel; for 
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TABLE V T 
CRITICAL VALUE OF K, AND CORRESPONDING VALUES OF K, c¢) 
Aluminium Mild loadi 
alloy weet chromium ype of loading Cc 
Steel st 
Critical value of K, 8°5 6 11 Direct stress and tl 
rotating bending cI 
Maximum value of K, cor- 6°5 4:25 9 Direct stress n 
responding to above K, p 
value 4 3:25 6°5 Rotating bending | 
Limiting value of K,=K,* 55 Le 8 Direct stress 
2 zs 5 Rotating bending 

b 
instance, with the nickel chromium steel in reversed direct stress the limiting K,=K;, Pp 
value is 8, which is more than twice the corresponding value for the mild steel. ) 

Hence it can be concluded that within the range of experimental conditions : — 
(a) The strength reduction factor K; equals the geometric factor K, up to 6 

some limiting value of K,. This limiting value of K, (=K;) varies with the 
material and the type of fatigue stressing employed (direct stress or i. 


rotating bending). 


(b) The value of K; reaches a maximum at some critical value of K;. 
This is again different for different materials, but the critical value of K; 
appears to be the same for direct stress and for rotating bending, although 
the corresponding maximum value of K;, is greater for reversed direct stress 
than for rotating bending. 


(c) Below the critical value of K,, non-propagating cracks could not be 
detected. However, for notches with a K, value above the critical value, 
non-propagating cracks were always formed in all materials and for both 
methods of testing. Although photographs in this paper are only given 
for cracks formed at stresses just below the fatigue limit, non-propagating 
cracks can be formed at much lower stresses, as mentioned in Refs. 12 
and 13. 


(d) Experimental points taken from other published work show general 
agreement with the present results. 


5. Discussion on the Relation Between K, and K, 


It has been shown that at a certain critical value of K, a non-propagating crack 
can form at the bottom of the notch. The shape of the curves in Figs. 25 and 26 
suggests that the relationship between K, and K;, to the left of the line marking the 
critical K, value is completely different from that existing to the right of the line. 


* Above this value K,< K, 
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The possibility that there may be two separate parts to this relationship has been 
overlooked by many investigators. 


It is suggested that this discontinuity is due to the formation of non-propagating 
cracks. Below the critical K, value, the relation between K, and K; depends on the 
stress necessary to initiate a fatigue crack; the nominal stress at the fatigue limit is 
that necessary to initiate the crack, the propagation being automatic. Above the 
critical value, a crack will form at a nominal stress below the fatigue limit; and the 
nominal stress at the fatigue limit of the specimen is that which is necessary to 
propagate the crack. 


It is therefore improbable that a universal relationship between K, and K;, exists. 


Although in Table V a limiting K,=K; value is given, it is not implied that 
below this value K; always equals K,. Size effect and invalidity of the maximum 
principal stress criterion of failure will tend to make K;, less than K,. The accuracy 
of the fatigue limit determinations due to testing only a few specimens at each K, 
value is insufficient to resolve these points in the present work. 


6. Conclusions 


The following conclusions apply only in the case of completely reversed loading, 
ie. Zero mean load. 


(i) Non-propagating cracks can form at the root of notched fatigue specimens 
when a critical value of the geometric stress concentration factor K; is 
exceeded This value varies for different materials. 


(ii) Because of the formation of these cracks it is suggested that there are 
two distinct parts to the relationship between K, and K;. Below the 
critical K, value, the nominal stress at the fatigue limit is that necessary 
to initiate the crack, propagation then being automatic; above the critical 
K, value, a crack will form at a nominal stress below the fatigue limit, the 
applied nominal stress at the fatigue limit being that necessary to 
propagate the crack. 


(iii) At the critical value of K, the strength reduction factor K; reaches a 
maximum value. While the critical value of K, is the same for both 
rotating bending and reversed direct stress, the corresponding values of 
K, may be different. 


(iv) For values of K, above the critical value, the value of K; does not increase 
and may tend to decrease with further increase in K,. 

(v) In practice, the increase in strength of specimens containing sharp notches 
which occurs when a non-propagating crack is present may not be 
realised, because the component will be subject to additional loadings, 
such as superimposed mean tensile loads, impact effects, and so on. 
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Some Exploratory Three-Dimensional 
Jet-Flap Experiments 


J. WILLIAMS and A. J. ALEXANDER 
(National Physical Laboratory, Teddington) 


Summary: This paper describes a first attempt to establish the order of magni- 

tude of finite aspect ratio effects for jet-flap wings. Lift, pitching moment, 

and drag results are presented from pressure-plotting measurements on a jet- 

flap model of aspect ratio slightly less than three. These experiments were 

complementary to two-dimensional tests carried out at the National Gas 
Turbine Establishment. 


1. Introduction 


The development of the jet-flap scheme at the National Gas Turbine Establish- 
ment has already been discussed by Davidson". The basic theoretical concepts 
and the early two-dimensional experiments have been described by Stratford’’’ and 
Dimmock™? respectively. In association with these researches, three-dimensional 
wind-tunnel tests were carried out at the N.P.L., to obtain quickly a first idea of the 
magnitude of finite aspect ratio effects for jet-flap wings. 


A small pressure-plotting model of rectangular plan form, with a 124 per cent 
thick elliptic section and a jet angle of 31-3°, was mounted with a plate at one end 
to give an effective aspect ratio slightly less than three. The lift, pitching moment 
and thrust forces were evaluated from measurements of the surface static pressures 
and the jet reaction. These tests were directly complementary to the N.G.T.E. 
two-dimensional tests”, in that the models were essentially of the same shape, size 
and general construction. 


The present paper, after outlining the experimental arrangement used in the 
N.P.L. three-dimensional tests, mainly describes the force results obtained on the 
wing. Details of the wing-section pressure distributions and of the general flow 
characteristics are described elsewhere“. 


NOTATION 
a, slope of two-dimensional lift-incidence curve 
A_ wing aspect ratio 


c wing chord 
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PRESSURE 
TURNTABLE, 
NS 
AIR FLOW 
AIR 
TO BROOMWADE RECEIVER 
AIR COMPRESSOR 
ik CONTROL (IN-INSIOE OAAMETER) 
{ 
‘FLEXIBLE HOSE 
Gam) (21N.INSIDE DIAMETER) 
Fig. 1. 
Arrangement of the jet flap model and external ducting. 

Ciy,Cmp»Cpp pressure lift, pitching moment and drag coefficients, (about 
mid-chord); derived by integration of pressure forces on wing 
surface 

Ci,Cm,Cr total lift, pitching moment and thrust coefficients, (about mid- 
chord); derived by adding direct jet reaction components and 
pressure forces on wing. (See Section 2.3) 
C; jet reaction coefficient=J/(4 p,U,7S) 


e 


aspect ratio efficiency factor 


J jet reaction at nozzle 
Po» Po» U, main stream pressure, density and speed 
p,, duct stagnation pressure 
R_ main stream Reynolds number based on wing chord 
s span of wing (excluding small ellipsoidal tip) 
S area of wing plan form (excluding small ellipsoidal tip) 
w width of jet nozzle 
a wing incidence relative to main stream 
6 jet exit angle relative to wing chord line. 


2. Experimental Method 


2.1. MODEL AND TEST SET-UP 


The model was of 8 in. chord and 12 in. span, with a small ellipsoidal tip fairing 


which added a further 0-25 in. to the span and 14 per cent to the wing area. It was 
mounted centrally in a near-circular tunnel of 7 ft. diameter so that tunnel inter- 
ference effects were negligible. The general arrangement of the model in the 
tunnel, and of the external ducting to the compressor unit, is illustrated in Fig. 1; 
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JET FLAP EXPERIMENTS 


a turn-table attached to the model root rotated in a fixed end-plate as shown. Mainly 
because of structural and accessibility considerations, the end-plate size was limited 
to one chord height vertically above the zero incidence chord line and two chords 
below, to 0-8 chord upstream of the wing leading edge and two chords downstream 
of the trailing edge (jet exit). According to Mangler’s end-plate theory’, the 
effective aspect ratio of the wing was then about 2°75, for the ratio (aerofoil span)/ 
(end-plate height) equal to 4. This effective aspect ratio is only slightly less than 
the value 3 which would have been obtained with an infinite end-plate. 


The chordwise pressure distribution was measured at each of the four spanwise 
stations, 0-2s, 0-5s, 0-8s and 0:95s from the wing root, by about 26 static pressure 
tappings—most closely spaced near the leading edge and the trailing edge. The 
blowing slot was located at the trailing edge, and had a mean width w=0-025 in. 
(w/c=0-:003), with a spanwise variation of less than 0-0025 in. 


At zero main stream speed, the jet emerged at an angle of 31-3° to the chord 
line for the whole range of jet total pressures pp measured in the model duct (up to 
two atmospheres absolute). The jet reaction J was determined by the N.G.T.E., in 
terms of pp, from thrust balance measurements at zero main stream speed; a 
correction was applied to allow for the surface pressure drag arising on the model 
as a result of the flow induced by the jet efflux. 


2.2. RANGE OF EXPERIMENTS 


Experiments covering the incidence range —5° to 20° were mostly at a wind 
speed of 100 ft./sec. (R=0-4 x 10°), when the available air supply permitted values 
of the jet reaction coefficient C; (=J/(4,,U,"S)) up to 0°5. Higher C,-values were 
obtained by reducing the wind speed to 50 ft./sec. (C; << 2-1). For a few tests, at 
zero incidence only, the speed was reduced still further to 30 ft./sec. giving C,-values 
up to 4-8. Transition wires were located on the upper and lower surfaces of the 
model at 0-2c behind the leading edge, as far forward as possible without causing 
interference at the closely spaced static holes in the wing nose. Except at very low 
C,-values, the results were not significantly altered by removal of the transition 
wires. 


2.3. DETERMINATION OF FORCE COEFFICIENTS 


The sectional “ pressure-force” coefficients were evaluated by chordwise 
integration of the surface static pressures at each spanwise station. Some simplify- 
ing assumptions were necessary close to the trailing edge, since there were no static 
holes inside the slot throat. The resulting error in lift is negligible, but this may 
not be so for the pitching moment and drag under all conditions. The wing 
pressure-force coefficients C1,, Cm, and Cp, were finally determined by spanwise 
integration of the sectional pressure forces. 


The total force on the wing was derived by adding the jet reaction at the nozzle 
to the pressure force on the wing arising from the main stream flow over the 
aerofoil surface. 
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Fig. 2. 


Spanwise distribution of pressure lift at zero incidence. Variation with C,. 


Thus, the overall lift, pitching moment, and thrust coefficients are written as 


CL sin (6 + a)+Crp 


The symbol 7 here denotes the perpendicular distance, from the point about which 
moments are taken, on to the extended centre line of the jet exit; for moments about 


C,=C, cos (6+ 2)—Cp» 


mid-chord, as in the present tests, n/c=4 sin 0. 
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SYMBOL U,FYAEC) 4 
° 100 
x 50 
4 30 
20 


PRESSU 


--~| JET REACTION 
LIFT | 


0-2 0-4 0-6 0-8 ye 12 4 1-6 20 
Fig. 3. 


Variation of total lift C, with C,} at zero incidence. 


3. Experimental Results 


In order to assess quickly the importance of finite aspect ratio effects on 
jet-flap wings, the experimental results on the present wing, of aspect ratio 2°75, are 
discussed in terms of simple aspect ratio considerations for conventional wings 
(without blowing) at incidence. 


3:1.. 


The spanwise distribution of the pressure-lift loading induced by blowing from 
the trailing edge, with the wing at zero incidence, is little different from that given 
by simple lifting-line theory for the plan form, or by experimental distributions for 
the wing at incidence without blowing (see Fig. 2). 


The total lift coefficient C;, at zero incidence is plotted in Fig. 3 against C,?. 
” At C,-values below unity, the experimental results satisfy approximately the 
relation 


)  —)} and are about 0°6 of the corresponding N.G.T.E. two-dimensional values. The 
| appropriate value of the aspect-ratio correction factor [1+(a,/(~Ae))]~' for lift on 
/ conventional wings is 0-58, assuming ideal elliptic loading (e=1) and a sectional 
h lift-curve slope a, of 2x. At higher C,-values and incidences (below the stall), the 
t tate of increase of C, with C,} tends to rise slightly. When C, > 2, the jet reaction 
) contribution to the total lift exceeds the pressure-lift contribution. 
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Fig. 4. 


Variation of total lift with incidence at constant C,. 


Lift-incidence curves for a range of C,-values up to 2:1 are shown in Fig. 4. The 
measured dC,/dz, at constant C,; and small incidences, increases from 0-055 /degree 
without blowing to about 0-1/degree when C;=2. The increase is roughly propor- 
tional to C;* and comprises contributions from both the jet reaction and pressure- 
force components. If a is small and measured in degrees, equation (1) with 
6=31-3° gives dC,/dz~ 0-015 C, + 
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Fig. 6. 
Variation in position of centre of total lift. 
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There is no significant loss in stalling incidence as C, is increased from zero, and at 
C,-values above unity there is even some improvement. However, this may be 
peculiar to the particular wing configuration of these tests. 


3.2. PITCHING MOMENTS 


The total pitching moment C,, about the mid-chord axis is plotted against the 
corresponding total lift C; in Fig. 5, for prescribed values of 2 and of C;. The 
mean slope (dC,,/dC,) of the curves for constant C, (a varied) is about 0-25 at 
low but non-zero values of C,*, so that the aerodynamic centre is then located close 
to the quarter-chord position. The slope decreases as C, increases, so that the 
aerodynamic centre tends to move further aft, by about 0-04c as C, is raised from 
0:2 to unity. 


As C; is increased with 2 constant, or « is decreased with C, constant, the 
nose-down value of C,, becomes steadily larger, so that the centre of total lift then 
moves rearward. Fig. 6 shows that this chordwise movement, away from the 
mid-chord position for z=0 and C, small, is quite appreciable. 


3.3. ToTAL THRUST AND PRESSURE DRAG 


Values for the total thrust coefficient C; on the wing can be determined, as in 
Fig. 7, by subtracting the measured pressure drag Cp, on the aerofoil surface from 
the corresponding component C, cos (?+ 2) of the jet reaction. This does not, of 
course, include the skin-friction drag. 


Now, for two-dimensional inviscid flow in the main stream and the jet, 
Cr=C,; thus, in addition to the direct jet reaction component C,; cos (6+), a 
recovery of thrust C, {l1—cos(@+2)} should theoretically be manifest in the 
pressure distribution over the aerofoil surface. The pressure drag measured in 
these three-dimensional experiments may therefore be regarded as comprising a 
section form drag, and an induced drag resulting from downwash effects at the wing 
less any thrust recovery which might arise. Because of the low aspect ratio and 
jet angle, the combined form and induced-drag contributions completely outweigh 
any thrust recovery (see Fig. 7). Moreover, with the present simple wing section 
and jet nozzle shapes, the recovery of thrust is probably small in practice. 


The measured pressure drag on the wing at zero incidence, with C,; <0°5, 
satisfies approximately the relation 


Cp, =0°013+0°14C,,”, . ‘ (3) 


while for C,-values up to 2, Cp, does not exceed 0-013 + 0-16 C,,” (see Fig. 8). An 
analysis of data for higher incidences up to the stall gives much the same result. 
The corresponding induced (pressure) drag C,,”/(*Ae) for a conventional wing 


*The unusual slope obtained for dCm/dC,, without blowing is ascribed to the unconventional 
trailing-edge shape and low Reynolds numbers of the tests. With a small amount of blowing, 
sufficient to clean up the flow near the trailing edge, the more usual value of 0°25 was found. 


28 The Aeronautical Quarterly 


| 
zz 


at 


JET FLAP EXPERIMENTS 


| | ‘ 
IDEAL THRUST RECOVERY a 
€,(-cos@) 
, 41 
1-25 
SYMBOL U,(FT, SEC.) 
° 100 ft 
x 50 
4 30 | PRESSURE ORAG 
Cop 
c 
T 97 
0-75 
A 
JET REACTION COMPONENT 
c,cos 
050 
iy 
0-25} 
0-25 0-50 0°75 1-00 1-25 1-50 75 2-50 
Fig. 7. 
Variation of total thrust C,, with C, at zero incidence. 
0-30; | xz 
\ 
0-25 = — 
SYMBOL U, (FT/SEC) \ 
° 100 \ C,-2-08 
x 30 
0-20 3 
“op 4 
Op tp 
0-15 — 
C,=0°49 L 
0-10) >— 
30-27 
0-05 
of 
025 0:50 075 +100 1-25 1-50 
Cis 
Fig. 8. 


February 1957 


Variation of pressure drag C,,, with C,,* at zero incidence. 
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producing the lift C,, is 0-116 C,,”, assuming elliptic loading (e=1). The difference 
0:013+0-02 C,,? from the value given by equation (3) could reasonably be 
accounted for as sectional form drag in the experiments. In any case, it seems that 
the pressure drag on a jet-flap wing is closer to that of a conventional wing 
producing the same pressure lift rather than the same total lift. 


4. Concluding Remarks 


The experimental results for this wing of aspect ratio 2°75 demonstrate the 
order of magnitude of finite aspect ratio effects on jet-flap wings and can usefully 
be reviewed in the light of simple aspect ratio correction factors for conventional 
wings at incidence (without blowing). Some supplementary experiments’, with 
the aspect ratio of the jet-flap wing increased to nearly 7, also lend support to the 
comparisons made in Section 3. However, the present tunnel tests were not 
designed to establish basic principles for a three-dimensional jet-flap theory, nor to 
resolve fundamental downwash considerations*. For this reason, and to provide 
more comprehensive data for project work, further experiments are to be made at 
the N.P.L. on a larger-scale model with variable aspect ratio and jet angle. The 
forces will be determined both by three-component balance measurements and by 
pressure-plotting. The three-dimensional nature of the flow over the wing and in 
the wake will also be explored. 


ACKNOWLEDGMENTS 


The work described in this paper was done in the Aerodynamics Division of the 
National Physical Laboratory, and is published on the recommendation of the 
Aeronautical Research Council and by permission of the Director of the Laboratory. 


REFERENCES 


1. Davipson, I. M. The Jet Flap. Journal of the Royal Aeronautical Society, January 1956. 


2. STRATFORD, B. S. Early Thoughts on the Jet Flap. The Aeronautical Quarterly, Vol. VII, 
pp. 45-59, February 1956. i 


3. Dimmock, N. A. An Experimental Introduction to the Jet Flap. National Gas Turbine 
Establishment Report R.175, July 1955. (To be published in the A.R.C, Current Paper 
Series.) 


4. WILLIAMS, J. and ALEXANDER, A. J. Three-Dimensional Wind-Tunnel Tests of a 30° Jet 
Flap Model. A.R.C. Current Paper 304, November 1955. 


5. MANGLER, W. Die Auftriebsverteilung am Tragfliigel mit Endscheiben. Luftfahrtforschung, 
Vol. 14, p. 564, 1937. 


*Some preliminary remarks on possible fundamental postulates are made in Ref. 4, but these 
are largely superseded by more recent theoretical studies at the Royal Aircraft Establishment, 
Farnborough. 


30 The Aeronautical Quarterly 


The Aiming Problem in Aerial Gunnery 


M. HANCOCK, A.R.CS., B.Sc. 


(Royal Aircraft Establishment) 


SumMMaRY: The aiming problem in aerial gunnery has formerly been analysed 
by approximate methods which allow results of sufficient accuracy for practical 
requirements to be obtained from a very simple mathematical treatment. 
There has therefore been no practical demand for any more elaborate analysis, 
although the development of a more exact treatment of the problem is not 
without interest from the academic point of view. In this paper a more 
complete analysis is developed in a form which is generally applicable to a 
broader class of aiming problems. The relevant parts of the standard external 
ballistic analysis for projectiles fired from guns are also included for the sake 
of completeness, and the earlier simplified analysis of the aiming problem is 
also presented to allow comparison of the results with those of the more 
complete analysis. Finally, a section on the gunner’s view of the aiming 
problem is included, to provide an explanation of certain difficulties which 
frequently arise through neglect of the important differences between the 
points of view of an observer moving with the gun and an observer at rest 
in the air. 


1. Introduction 


It is well known that in aerial gunnery the gun must generally be pointed 
slightly away from the line joining it directly to the target to obtain a hit, the angle 
by which it is turned away being called the “lead angle”. The exact calculation 
of the lead angle in the great variety of circumstances which can arise in air combat 
is a matter of considerable complexity, and no complete analytical treatment of the 
problem has been published previously. From the point of view of those concerned 
with the design or use of aircraft gunsighting systems this is not a serious matter, 
as analytical results of adequate accuracy for all practical requirements can readily 
be obtained by treating the problem in a number of separate parts which are 
regarded as independent, although they are in fact interdependent to some extent. 


It is of some interest, however, to examine the problem analytically as a single 
entity instead of in separate supposedly independent parts, and this paper presents 
the results of such an analysis. In this new analysis simplifying assumptions 
equivalent to supposing the action to take place over a flat earth in an atmosphere 
of uniform density free from velocity gradients are still retained, but these condi- 
tions are so nearly satisfied in actual combat that their introduction introduces 
errors of extremely small significance. 
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speed of the gunner’s aircraft along G,V at the moment the shot is 
fired 


direction cosines of G,V 


Vy muzzle velocity of the projectile 
R_ target range G,T, at the moment the shot is fired 
R rate of change of R 
(lz, Mr, Mz) direction cosines of the sight line 
6 the angle T,G,V between the line of sight and the gunner’s line of 
motion 
» angular velocity of the line of sight. Orthogonal components of w 


(ly, My, Ny) 


%o 


Pis Par 


denoted by ,, w,, and w,’, w,’, are also required 
direction cosines of the gun line G,G, 


the angle G,G,V between the gun line G,G, and the gunner’s line 
of motion G,V 


the lead angle 7,G,G,. Orthogonal components of z denoted by 
a,, % and a,’, a,’ are also required 


initial velocity of the shot along G,N 

time of flight of the projectile from G, to T, 

the distances G,N, NM and MT, respectively 

the projections of 7,’T, on OX, OY and OZ respectively 
magnitude of yaw angle 


orientation of the plane of yaw relative to an arbitrary reference 
plane 


initial values of 6 and @ 
a constant defined in Section 4.3 


functions defined in Section 4.3 


i= /(-1). 


2. General Description of the Analysis 


The analytical treatment of the aiming problem falls into two parts which are 


examined separately initially, although the results of the two parts must ultimately 
be combined. These two parts are (i) the Kinematic part, in which the geometrical 
relationships are examined, and (ii) the Ballistic part, in which the dynamic aspects 
of the projectile motion are examined. The new analysis presented in this paper is 
confined to the kinematic part of the problem, and this is presented first. The 
ballistic part of the problem has already been fully discussed in several places":*-*:* 


*Some duplication of symbols occurs between the notation used in the discussion of the 
kinematic aspects of the problem and that used in the discussion of the ballistic aspects of 
the problem, but no danger of confusion is introduced by this. 
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and the results, in so far as they are relevant to the present paper, are quoted with a 
minimum of discussion. The simplified analysis referred to in Section 1 is then 
presented for comparison with the new analysis, and finally a short description of 
the aiming problem as the gunner sees it is added. 


3. The Kinematic Aspects of the Problem 


3.1. THE GENERAL PROBLEM 

The kinematic analysis is completely general, and is applicable to the release of 
projectiles of any type, including rockets and bombs from aircraft in flight, so long 
as the trajectory elements of the projectile can be expressed in the form required. 
The analysis is referred to the general air combat situation illustrated isometrically 
in Fig. 1, where GG,G represents the path of the gunner’s aircraft and T,TT, that 
of the target. A correctly aimed shot fired by the gunner from G, when the target 
is at T, hits the target at T,, the direction of the gun line required to achieve this 
being G,G,. The gunner’s line of sight to the target at the moment the shot is fired 
is G,T,, and his direction of motion at the same instant is G,V. The lead angle is 
thus T,G,G,, and the initial direction of motion of the shot through the air resulting 
from its muzzle velocity along G,G, and the gunner’s own velocity along G,V is 
along G,N. The trajectory of the projectile is G,ST, which lies close to, but 
generally not exactly in, the vertical plane through G,N. The displacement of the 
shot between G, and 7, can be resolved into the three rectilinear components, G,N 
along its original direction of motion, NM ' at right angles to the vertical plane 
through G,N, and MT, vertically. The direction of motion of the target at the 
moment the shot is fired along T,7,’, tangential to 7,77, at T,, while the distance 
T,T,’ represents the target displacement during the time of flight of the shot due 
to the velocity of the target at T,. The distance T,’T, represents the component of 
total target displacement due to changes in target velocity between T, and 7). 


The analysis is referred to non-rotating three-dimensional right-handed 
co-ordinate axes with G, at the origin, OZ vertical, and OX and OY horizontal. 


3.2. THE KINEMATIC ANALYSIS. FIRST PART 


The kinematic analysis is carried out in two separate parts. In the first part 
three auxiliary construction lines are required in Fig. 1. The first of these, shown 
as G,P. is located in the plane T,G,V in such a position that the plane G,G,P is 


Fig. 1. 
The general air combat situation. 
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perpendicular to T7,G,V. The lead angle T,G,G, can thus be measured in terms of 
the two orthogonal components 7,G,P and G,G,P denoted by «, and 2, respectively. 
The second auxiliary line N, is normal to T,G,V and the third N, is normal to N, 
and to G,T,. These two normals are not shown in Fig. 1. They have direction 
cosines (I,, m,,n,) and (I,, m,, respectively, defined as follows :— 


1, =(meNr — Mang) cosec 


1, =(lg—Ip cos 9) cosec 6 
m,=(mg— Mp COS 9) cosec ‘ (2) 


nN, =(Ng— Np Cos 9) cosec 


The analysis is carried out by examining the conditions which must be satisfied 
for the shot and the target to arrive simultaneously at T,. The algebraic manipula- 
tions involved in this are of some complexity, and have already been published else- 
where”. The results take the simplest form for the special case where the transverse 
component d, of the trajectory is small enough to be neglected, a condition which 
is very nearly satisfied in most ordinary circumstances. For this case the components 
of lead angle are given by 


,t—d,) sin 0+ 


VV ,t—d,) cos + mpdy + + d;)} 
Vo 
and sin a, = ——— [Rwo.,t + 1,6x+m,dy +n,(6z+ d,)] ‘ (4) 
d Vy 


where , and », are the component angular velocities of the sight line about N, 
and N, respectively, and the quantities V,, t, d, and d, are inter-related by the two 
equations 


V2(V ot d,)?+2V t—d,)[(R + +Rw,tsin 6+ lgdx + medy + ng(8z+ d;)]+ 
+ Rt)? + (Rot)? + (Ro,t)? + 8x? + 8y? + (8z+d,)? + 


+2(R+ R1){Ipdx + mpdy + np(8z + d,)} + 2Ro,t{l,5x + + n,(6z+d,)} + 


and V 2d, -2V.Vol(R + Ri) cos 64+ Vet + Ro,t sin 6 + + + ng(dz + d;)] + 


and by two further equations derived from the ballistic analysis. 
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For the more general case where d, is not negligible, an approximate solution 
can be obtained by adding the quantities 


Vulz)? to dx 
and dV gla t+ Vulz)/((Vala + V ule)? +(Vama+Vumr)*}? to dy 


in equations (3), (4), (5) and (6), and by introducing a further relationship to define 
d, from ballistic considerations. 


An exact solution for the case where d, is not negligible is given by 


. _, [1 -sin?a,]}* _,B 
and 
(X, vX.){(manr— MpNg :) +V (nelr ] 
“Vad, (1+?) + — | cosec 
(8) 
where 


=[(R+ + Velet + cos 6) + Mpng)} cosec 
X,=[(R+Ri)mp + Vamet + Rt{o,(me— mz cos cosec 6 + 8y] 
X,=((R+ Rt)np + Vgngt+ Rt{o,(ng— Np Cos 6) cosec 6+5z+d,] 


{& cos (ma — mp cos 


+ X,(ng — Np COS - 


— Vg sin? cosec 6 


_ 1 FV (X,+vX,) VoX ] 


and V,, t, d,, d, and d, are inter-related by the equations 


+ 0X2)? + (1+ 2d, VVal(X, + (la + vig) + (1 + v*) + 
. (10) 


and V,7-2 (1 


+moX, | +e? - Vu2)=0 an 


and by three further equations from the ballistic analysis. 
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3.3. DISCUSSION OF THE FIRST PART OF THE KINEMATIC ANALYSIS 


The kinematic analysis outlined in Section 3.2 gives equations which, when 
combined with the results of the ballistic analysis, allow the required direction of 
the gun line to be determined from the primary independent variables to as close 
an accuracy as is consistent with the inherent uncertainties associated with the 4 
terms. Apart, therefore, from the extreme complexity, the equations defining the 
components of lead angle could be used for the design of a gun sighting system 
having the line of sight directly controlled by the gunner and the gun direction 
controlled automatically in accordance with the equations. Sighting systems of this 
type, known as “ director systems,” are in common use for the control of heavy 
guns in ship and land based installations, but their complexity and cost are too 
great to make them generally suitable for use in aircraft installations. In aircraft 
a simpler form of mechanism known as a lead-computing sight is almost universally 
employed. In this the gunner exercises direct control, not over the sight line but 
over the gun line, the sight line being linked to this and following its movements 
through a quasi-elastic coupling in which the computations of lead angle take place. 
In a typical lead-computing sight, the sight line is defined by purely optical means 
and has no tangible mechanical existence, with the result that the magnitude of 
the angle @ and the plane containing it cannot readily be determined, and there is 
therefore an important gap in the information required in the equations of 
Section 3.2. 


The lead-computing sight does, however, allow information to be obtained 
about the angular velocity of the sight line, and it allows the angle ¥ and the plane 
containing it to be determined, besides providing a measure of the lead angle. To 
obtain results which can be applied directly to the design problems of such a sight, 
the kinematic analysis must therefore be approached in a manner which treats y 
and the plane of ¥ as independent variables, while @ and the plane of 6 are treated 
as unknown quantities. The analysis must also make use of the known components 
of angular velocity of the sight line, and the results will then show the components 
of lead angle referred to the Y-plane. 


The relationship between the planes and angles associated with this argument 
is shown in Fig. 2 which shows the gunner’s line of flight G,V, the line of sight to 
the target G,T, and the gun line G,G,. The line G,P is located in the plane 
T,G,V (the plane of 4) in such a position that the lead angle components T,,G,P (z,) 
and G,G,,P (z,) are in mutually perpendicular planes. In the analysis of Section 
3.2, G,T, and G,V are given, the angle 9 between them and the plane of @ being 
known, while the required position of G,G, is determined from equations which 
define the angles 2, and a, in terms of the angular velocity of G,T, and other 
measurable quantities. In the lead-computing sight, however, the angle 6 and the 
plane containing it cannot be found, although the angle ¥ between the gun line and 
flight line and the plane containing ¥ can be found, while the angular velocity of 
the sight line can also be determined. In this situation the lead angle is defined 
by the orthogonal components z,’, z.’, located in and at right angles to the plane 
of ¥. For this purpose an additional construction line G,Q in the ¥-plane must be 
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Fig. 2. 
Relationship between details of first and second kinematic analyses. 


introduced to make the component lead angles G,G,Q (z,’) and T,G,Q (z,’) mutually 
perpendicular. The kinematic analysis associated with the determination of 2,’ 
and 2,’ is outlined in the next section. 


3.4 THE KINEMATIC ANALYSIS. SECOND PART 


In the second part of the kinematic analysis, the angular velocity of the sight 
line is measured in terms of the components ©,’, which produces motion of the 
sight line at right angles to the plane 7,G,Q, and »,’, which produces motion of 
the sight line in the plane T7,G,Q. Three auxiliary construction lines are required: 
N’ perpendicular to the plane G,G,V; N,’ about which the component angular 
velocity »,’ occurs; and N,’ about which the component angular velocity ,’ occurs. 
The analysis is carried out as in the first part by examining the necessary conditions 
for the shot and target to arrive simultaneously at T,. The results are as follows: — 


Rwo,’t 


7 -1 
(R+Rt) +(R+Rty} 
where 
Y= - d,)Vq cos Y—d,Vy] + [ludx + mydy + ny(5z + + 
0 


(Valu + Vela)? +(V + Veg)? 
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Y= -d,)Vq sin 
Vo 


+[(lg— ly cos + (mg my cos + (Ng Ny Cos (6z + d,)] cosec 


d(lqmy — lyme) (Vu + Va Cos cosec 


(Valu + Valo)? +(Viurnu + Varna)? (15) 
Z =[(myne Meny)x + (Mula + (uma — (6z + cosec + 
d[Vulna nu Cos ¥) + Volta cos - ny)] cosec 


(Valu + Vale)? + + Vax)? 


and the relationships corresponding to equations (5) and (6) of the first analysis are 
provided by 


and (R+Rt)?+(Ro,'t? . (18) 


It is of interest to note that the results of this analysis are not complicated by 
the presence of the d, term to the same extent as in the results of the first analysis. 


3.5 APPROXIMATIONS DERIVED FROM THE KINEMATIC ANALYSIS 


The results obtaiied in Sections 3.2 and 3.4 are too complex for practical 
aj plications, and it is of interest to examine the approximate results, which can be 
obtained by expressing these results as series in ascending powers of small terms. 
The following results are obtained from the first analysis : — 


a= (t- sin 8+ ,t+ z [ 18x + m,dy + n,(8z + | +higherterms (19) 
and @,= + 1,6x+m,dy+n,(6z+d,) | +higher terms. (20) 
Vut Vuty 


where f, represents the ratio d,/V,, i.e. the time the projectile would take to travel 
the distance d, if it travelled without retardation. 


The corresponding approximations from the second analysis are 


,_ Rot Va(t ) 
Vute + 1 ) sin + higher terms (1) 


The results given in equations (19) to (22) will be found to bear a close resem- 
blance to the results of the simplified analysis presented in Section 5. They illustrate 
by their form that the components of the lead angle can be treated, to a close 
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approximation, as the vector sums of independent parts, depending on the angular 
velocity of the line of sight and the retardation of the projectile due to air resistance. 
The gravity components of lead angle are represented by the terms in d,. 


4, The Ballistic Aspects of the Problem 
4.1 THE GENERAL PROBLEM 


The ballistic part of the aiming problem is concerned with the factors which 
influence the motion of the projectile along the trajectory, and with the development 
of expressions showing the relationship between the time of flight 7, the trajectory 
projections d,, d,, d, and the initial conditions under which the projectile is fired. 
This subject is already covered by a considerable literature"-*-*-*) and it is only 
necessary here to quote those results which are of most significance in aerial gunnery. 


The whole pattern of the ballistic problem and its solution is set by the fact 
that the resistance of the air to projectile motion cannot be represented exactly by 
any mathematical expression, and therefore no exact solution of the equations of 
motion of the projectile can be obtained. A number of approximate procedures 
have been described, however, by which results of adequate accuracy can be 
obtained. The simplest of these which is suitable for all ordinary problems in aerial 
gunnery is based on an expression for air resistance of the form (Ap/c) v'*, where 
A is a constant which can be adjusted to give the best approximate fit over the 
velocity range considered; p is the air density, c is a constant which takes account 
of the projectile size and shape, and v is the velocity of the projectile through the air. 


In the most general trajectory analysis, account must be taken of the fact that 
the projectile may be yawed, i.e. its axis of figure may not point along its direction 
of motion, in which case the aerodynamic drag will be greater than it would be 
without yaw, and its direction may not be exactly in line with the direction of 
motion. The trajectory analysis with yaw is considerably more complex than 
without yaw, and it is preferable to consider the results for the simple case without 
yaw before examining the complications introduced by yaw. 


4.2. THe TRAJECTORY WITHOUT YAW 


When yaw is absent the transverse displacement d, is zero, and the relation- 
ships between ¢,, d,, and d, derived from the simple drag law quoted in 
Section 4.1. are 


V ot 
(1+J?) (23) 
1+4Jt 
2 
A pv V, 


and M denotes the mass of the projectile. 
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GUNNER'S 

AIRCRAFT 

VELOCITY Vg 

INITIAL PROJECTILE 
VELOCITY Vv, 


INITIAL ANGLE 
OF YAW 


VELOCITY 


PROJECTILE INITIALLY POINTS ALONG Vu BUT 
MOVES ALONG Yo 


Fig. 3. 
Initial yaw produced by firing into a cross wind. 


4.3. DESCRIPTION OF THE YAWING MOTION 


In general there are a number of causes by which yaw can be introduced into 
projectile motion. The most important of these in aerial gunnery is that produced 
by firing into a cross wind when the initial direction of relative motion between the 
air and the projectile is at an angle to the axis of the projectile, as shown in Fig. 3. 


The yawing motion produced by firing into a cross wind is described by the 
equation 


1@-9) 


2 


where 8 is the magnitude of the yaw angle 


¢ is the orientation of the plane of yaw relative to an arbitrary 
reference plane 


are the initial values of 6 and 


o is a constant associated with the dynamic stability of the 
projectile, having a value generally in the range 0:5 =o <1 


P:,P29:»42 are functions of the aerodynamic constants, and of the spin 
of the projectile, and of time. They increase monotonically, 
but not necessarily linearly with time from zero at the origin 
of the trajectory. 


Equation (26) represents a motion which is the sum of two damped rotating 
vectors having in general different initial magnitudes, different rates of rotation, and 
different damping factors. An illustration of the characteristic rosette pattern 
obtained by plotting 5 against ¢ in polar co-ordinates is shown in Fig. 4 for a typical 
set of circumstances. In this example the first two complete rotations of ¢ about 
the origin would occupy about the first 0-05 second from the start of the trajectory, 
and the magnitude of the yaw would be reduced to one tenth of its initial value 
after about 0:25 second. 
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= SAAN 


PLOT OF YAW (6) AGAINST ORIENTATION OF YAW (9) 


Fig. 4. 
Typical example of yawing motion. 


4.4. THE EFFECTS OF YAW ON THE TRAJECTORY 


From the practical point of view, the importance of the yawing motion which 
occurs when a projectile is fired into a cross wind is not in the yawing motion itself, 
but in the effects of the yaw on the forward and transverse motion. These effects 
may be considered under two separate headings which are discussed in decreasing 
order of practical importance. 


4.4.1. The Effects of Yaw on the Forward Motion 


When the projectile is yawed, the projected area it presents in the direction of 
motion is greater than in the unyawed condition and the aerodynamic drag is there- 
fore increased. Investigations have shown that, up to yaw angles of about 16°, the 
increase of drag due to yaw is closely proportional to the square of the angle of 
yaw, so that the effect of yaw on forward motion can be taken into account by 
introducing a factor of the form (1+?) into the equations of motion before 
integration. 


The solution of the trajectory equations when drag due to yaw is introduced 
is considerably complicated, but some simplification can be introduced by examining 
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the square of the yaw amplitude given by equation (26). This is 


The periodic term in this equation has an amplitude which is comparable with the 
non-periodic terms, but the variations in drag which it produces fluctuate so rapidly 
that its effect is of no significance in practical trajectory calculations. It is sufficient 
therefore for practical purposes to calculate the mean square of the yaw over the 
periodic time of the oscillatory term in equation (27) and to use this mean square 
value in the trajectory calculation. 


The resulting expression for the effective mean square yaw at time t, denoted 


by 5? (1) is 
1\? -2,+4) ( 1\? -2 


Even with this simplification, direct integration of the trajectory equations is 
not possible, and recourse must be had to numerical methods of integration. 


4.4.2. The Effects of Yaw on the Transverse Motion 


In addition to the increase in drag produced by yaw, another effect, generally 
of small practical importance, follows from the component of aerodynamic force 
acting at right angles to the direction of forward motion. This transverse force 
analogous to the lift on an aerofoil produces a transverse acceleration of the 
projectile, following the rotations of the yaw plane, with the result that the path of 
the centre of gravity is changed from the nearly straight line it would follow in the 
absence of yaw to a helical form lying close to this line. For analytical purposes 
the resulting departure from the no-yaw path is considered in two separate parts, 
the first being a rosette form of spiral motion similar in general appearance to the 
yaw rosette shown in Fig. 4, and the second a steady deviation from the 
no-yaw path. 


The rosette component of the motion is described by the equation 


where r and W denote the radius and orientation of the displacement in the spiral 
from the corresponding position in the equivalent plane trajectory. The similarity 
of this motion to the rosette yawing motion is apparent from a comparison between 
equations (26) and (29), but the first term in the square bracket of equation (29) is 
negligible in comparison with the second term, with the result that the motion takes 
the form of a simple damped spiral. In typical cases the initial radius of the spiral 
is only a small fraction of an inch so that this component of the motion is of no 
practical significance. 


42 The Aeronautical Quarterly 


ref 


pr 
of 


5, 
fo 
th 
a 
of 
fli 
be 
cc 
ps 
S 
Ol 
a 
i 
| 
tl 


AERIAL GUNNERY 


The second component of the transverse motion is a progressive displacement 
represented by the equation 


where H, and Z, denote the Cartesian components of the displacement. In most 
practical cases this transverse motion represents a departure of only a small fraction 
of a degree from the equivalent trajectory without yaw. 


5. The Simplified Analysis of the Aiming Problem 


In the simplified analysis of the aiming problem mentioned in Section 1 the 
following approximations are made. 


Firstly, it is assumed that the effects of gravity can be treated separately from 
the other factors which influence the lead angle. In making this simplification it is 
assumed that, when the position of the projectile has been calculated at any given 
time as though gravity were absent, its actual position will be at a vertical distance 
of 4gt? below that point, where g is the acceleration of gravity and ¢ the time of 
flight to the point in question. If air resistance were absent this assumption would 
be strictly accurate, and it is in fact sufficiently accurate for most ordinary 
conditions. 


Secondly, it is assumed that the target flies at constant speed along a straight 
path. Under ordinary conditions it is not possible for an aircraft to change speed 
sufficiently during the time of flight of the projectile to have any appreciable effect 
on the lead angle calculated in this way, although the acceleration which can be 
applied across the line of flight can curve its path sufficiently to have a serious 
influence on the lead angle at longer ranges. The effects of such curvature can be 
calculated with ample accuracy by first treating the target motion as being at 
constant speed in a straight line and then finding its displacement from the point it 
would have reached under these conditions by virtue of its acceleration. 


Thirdly, it is assumed that in the absence of gravity the projectile would travel 
through the air in a straight line whose direction is defined by the resultant of its 
muzzle velocity relative to the gun and the gun’s velocity through the air. This 
assumption is very nearly true, but it neglects the fact that the projectile moves 
through the air along a spiral path of small and rapidly decreasing amplitude, the 
axis of which is slightly inclined to the assumed straight path. 


The analysis is carried out by reference to Fig. 5 where G,GG, represents the 
path of the gun and 7,7, that of the target. The shot is fired when the gun is at 
G and the target at T,. The velocities are denoted by Vq for the gun, V» for the 
muzzle velocity of the shot and V; for the target. The initial velocity of the shot 
resulting from Vy and Vg is denoted by V, and the average velocity from G to T, 
by V,. The ratio of V, to V, is denoted by p and the time of flight of the shot 
from G to T, by t. The line of sight GT, makes angles @ and with G,GG, and 
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G, G G, 


Fig. 5. 
Construction for the simplified analysis of the aiming problem. 


T,T, respectively and GQ represents the gun line. The lead angle 7,GQ is 
denoted by z. Initially it is necessary to restrict consideration to the special case 
where G,GG, and 7,7, are coplanar, but it is shown later how the treatment can 
be extended to meet the more general case where this condition is not satisfied. 


To determine the lead angle, the point Q is found by drawing through 7, 
parallel to G,GG, to meet the gun line in Q. It is also necessary to draw T,N and 
MQL perpendicular to GT, and T,L parallel to GT,. The time which the shot 
would take to reach JT, in the absence of retardation, i.e. the vacuum time, is 1). 
The following relationships hold : — 


where MQ=GQsina 
NT,=T,T, sin 
QL=T,Q sin QT,L. 


Since and GT, =tV,=1,.Vo, 


GQ=1t,Vy. 


Also 7,7,=tV>. 
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Since GO 
GQ=1,Vo. 
Vu 
Also QT,L=6. 


Thus equation (31) is equivalent to 
t.Vy sin sin ¢—t,Vq sin 


which, on replacing sin a by 2 and dividing throughout by t,V., becomes 


a=p sin 6 ‘ (32) 


The angular velocity of the line of sight to the target at the moment the shot 
is fired is given by » R=V;sin¢-V, sin, and this, in combination with (32), 
leads to 


a= Dp sin 6 . 

In the more general case, where G,GG, and 7,7, are not coplanar, the result 

of equation (33) still applies, provided that z, w and Vq are regarded as vector 

quantities and are compounded accordingly. Alternatively, equation (33) may be 

regarded as applying only to components of w and a lying in the plane defined by 

the gunner’s path G,GG, and the sight line GT,._ A second equation from which 

the V, term is absent is then required to relate the components of 2 and » in a plane 
at right angles to this. 


In practical cases of aerial gunnery a further component of lead angle must 
also be introduced to allow for the gravity drop of the shot, which has been 
neglected in the foregoing treatment. This component lies in a vertical plane and 
must be calculated to compensate for the vertical distance which the shot falls under 
the influence of gravity during its time of flight to the target. 


The similarity between the result expressed by equation (33) and the 
approximations derived from the full kinematic analysis in Section 3.5 will be 
immediately apparent. 


6 The Aiming Problem as the Gunner Sees It 
6.1. THE GUNNER’S VIEWPOINT 


In the foregoing discussion the aiming problem has been examined from the 
point of view of an observer at rest in the air through which the gunner and target 
are moving. This treatment provides all that is strictly necessary to obtain a 
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solution of the problem, although it is useful in conclusion to examine the situation 
from the gunner’s point of view as well, because the gunner always feels himself to 
be at rest and may find some difficulty in reconciling what he sees with the results 
of the foregoing treatment, unless the matter is correctly explained to him. 


When tracer ammunition is used, the gunner can often see his bullets as they 
travel towards the target. For the first few yards they appear as streaks of light, 
but for the remainder of their flight he sees them as spots of light which generally 
seem to be moving on strongly curved paths. To the gunner’who knows that, apart 
from the effect of gravity, they travel through the air in straight lines, their appear- 
ance may be very surprising, because the observed curvature is often much greater 
than that due to gravity and may be in any direction. 


The apparent curvature of the bullet paths as seen by the gunner arises partly 
from the effect of gravity, but is also influenced strongly by the deceleration due to 
air resistance acting in combination with the gunner’s motion through the air, and 
by any angular motion of the sight line which may be present. 


6.2. CURVATURE DUE TO GRAVITY 


The influence of gravity is to curve the flight path of the bullet downwards 
slightly below the straight path which it would otherwise follow. This effect can 
be observed alone when a gun is fixed on the ground. The curvature which is 
produced is much smaller than that due to the other two causes. 


6.3. APPARENT CURVATURE DUE TO DECELERATION OF THE PROJECTILE 


To understand the second cause of apparent curvature it is useful to consider 
the special situation in which the gunner and target are travelling on parallel 
courses at the same speed. In the absence of air resistance the target would be hit 
if the gun pointed directly towards it, and the bullet would be seen by the gunner 
as moving directly towards the target. This ideal situation is shown in Fig. 6 
where G,G.G, represents the gunner’s path through the air, and 7,7,7, represents 
the target’s path. 


Ts Te Ti 


G3 Ge 


Fig. 6. 
Apparent curvature due to projectile deceleration. 
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The gun is fired at G,, while pointing towards the target at T,, and the shot travels 
along G,T, to meet the target at T,. The point G, represents the position of the 
gunner at any arbitrarily chosen instant during the interval between the firing and 
impact of the shot, and T, represents the position of the target at the same instant. 
The postion of the gunner at the moment of impact is represented by G;. In the 
absence of air resistance and gravity, the shot travels with constant velocity along 
G,T, and thus remains at the point S, where G,T, intersects G,T,, throughout its 
flight. The direction in which the gunner sees the shot thus remains constant along 
the gun line at the instant of firing. 


When air resistance is present the shot loses velocity as it travels along G,T;, 
and will thus only have reached some position such as S’ when the gunner and 
target are at the ends of the line G,T,. The gunner’s line of sight to the shot will 
therefore lag behind the direction of the gun barrel by the angle S’G,S, and this 
angle will increase as the shot travels away from the gunner, because the distance 
SS’ by which the shot lags behind its “no-resistance” position increases more 
rapidly than the distance G,S. The apparent path of the shot as seen by the gunner 
thus takes the form of a curve represented by G,S,’S,’ which trails behind the gun 
line, and the gun must therefore be aimed ahead of the target by the angle which 
this curve subtends for a hit to be obtained. This angle is known as the trail angle, 
and the distance T,S,’ is known as the trail distance. 


(2) (b) 
G G 
Fig. 7. 
Apparent curvature due to sight line rotation. 
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6.4. APPARENT CURVATURE DUE TO SIGHT LINE ROTATION 


The third cause of apparent curvature, rotation of the sight line, is most easily 
understood by considering what would happen in the absence of gravity and air 
resistance with the gun at rest and the sight line rotating. 


Fig. 7 (a) shows what actually occurs. The shot is fired from G when the 
target is at T, and travels along the straight line GS,S,. The shot positions S, and 
S, correspond to the same instants as the target positions T, and T,, and the shot is 
supposed, for the sake of the present argument, to pass as close as possible to the 
target without actually hitting it. When the shot first leaves the gun, the angle 
between the line of sight to the target and that to the shot, as seen by the gunner is 
T,GS.,, but when it reaches S, this angle is reduced to T,GS, and, as it passes the 
target, the angle is reduced substantially to zero. Afterwards, when the shot has 
passed the target and reached S,, the angle has changed sign and become S.G7T;. 


This situation as it appears to the gunner is shown in Fig. 7 (b). With his 
attention fixed on the target, the path of the bullet appears to him to curve in 
towards the apparently stationary line of sight to the target, and eventually to cross 
it. The angle through which the gunner’s line of sight must be turned to follow 
this curvature is given by the product of the angular velocity and the time of flight. 
In the general case it arises from the relative speed between the gunner and target, 
and it is therefore known as the relative speed component of the lead angle. 


6.5. THE GENERAL CASE 


In the general case the three causes of curvature will act simultaneously to 
produce the resulting curvature seen by the gunner, and the total lead angle will 
therefore be composed of a trail component in the plane defined by the sight line 
and the gunner’s line of motion, a relative speed component in the plane defined by 
the relative motion between the target and gunner, and a gravity component in the 
vertical plane. 


REFERENCES 


1. Cranz, C. and Becker, K. Handbook of Ballistics. H.M.S.O. London (Translated from 
the Second German Edition), 1921. 


N 


Fow ter, R. H., GALLop, E. G., Lock, C. N. H. and RICHMOND, H. W. The Aerodynamics 
of a Spinning Shell. Transactions of the Royal Society A. Vol. 221, pp. 295-387, 1920, 
and Vol. 222, pp. 227-247, 1922. 


3. Unpublished National Physical Laboratory reports. 
4. Unpublished Ordnance Board report. 
5. Unpublished Ministry of Supply report. 


48 The Aeronautical Quarterly 


A Finite Series Solution for Grillages under 
Normal Loading 


D. F. HOLMAN 


(Naval Construction Research Establishment) 


Summary: A method has been given by Klitchieff@: ?) for analysis of grillages 
of intersecting beams under lateral loading. Klitchieff’s solution was in the 
form of infinite Fourier series and the expressions for bending moments 
converged approximately as = n-*, making numerical computation difficult. 
This paper shows that, for grillages in which one set of beams are evenly 
spaced, Klitchieff’s solution can be adapted to give finite series for deflections 
and bending moments at the intersections of the beams. The number of terms 
in the series is not greater than the number of beams in the equally spaced set 
(referred to in this paper as set “ A’). A numerical example is given and an 
accurate solution is obtained with appreciably less computation than other 
known methods would entail. 


1. Introduction 


The two orthogonal sets of beams may be called sets “A” and “B.” In the 
procedure given by Klitchieff the deflection of the s‘* beam of set B is expressed as 
a Fourier series with unknown coefficients C,;, (s=1, 2,...... , q). Bending 
moments in the same beam can then be expressed in two ways, firstly in terms of the 
second derivative of deflection, and secondly in terms of the forces 
on the beam. By beam theory for the r* beam of set A, we can find the 
unknown reactions R,, in this second expression in terms of deflections at the 
intersections of this beam with those of set B. These deflections can then be 
replaced by their Fourier expansions. Equating the two expressions for bending 
moment and applying a finite Fourier transformation, a set of algebraic equations 
involving the q infinite sets of unknown coefficients is obtained; and, in cases where 
the beams of a set are uniform, equal and evenly spaced, by taking advantage of the 
periodic properties of these coefficients the equations can be reduced to a finite set 
with g unknowns. (If there is symmetry in the A-direction the number of 
unknowns is $q or $(q+1)). When these have been solved, the deflections and 
bending moments may be found as the sums of infinite series. However, the bend- 
ing moment expansion is of the form =n~*w (n), where w (n) is a trigonometrical 
function, so that computation presents considerable difficulties. The following 
adaptation of the method is intended to overcome these. 
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SETA 
(p BEAMS OF LENGTH b) 
rel / 


- ---------- 

SET B. 
LENGTH a) 


| 
| | 


(a) PLAN SHOWING NUMBERING OF BEAMS 
AND INTERSECTIONS 


Rs 
W, 


(b) ELEVATION OF BEAMS 


Fig. 1. 
The general grillage problem. 


NOTATION 


A,B _ the two sets of intersecting beams 
p,q numbers of beams in sets A, B, respectively 
w deflection 
F _ shear force 
M _ bending moment 
R,, reaction between the r beam of set A and s* beam of set B 
h= p+l 
a,b lengths of B-beams and of A-beams respectively 
I,,1, | moments of inertia of A-beams and B-beams respectively 
c,d spans between consecutive A-beams and B-beams respectively 
x distance from end of B-beams 
x, distance of r A-beam from end of B-beams 
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Coefficients of the deflection expansion 
Lue (mr)*Cns 
V = | m+2jh |, where j is any integer. 
Suffixes 
r,s refer to the r™ beam of set A and the s beam of set B. 


2. Conversion of Solution to Finite Series 


2.1. First CASE: SIMPLY-SUPPORTED ENDS 


The deflection of the s‘* beam of B can be represented by 
oo 
and the bending moment is then 


M,=EI, ax? = El, > a e (2) 
R,, is defined as the downward reaction on the s** beam at its intersection with 


the r beam of set A, so that, taking moments (if there is no other load on the 
B-beam), 


Res (1 + Ry (x—x,) U (x—x,), (3) 
where U@=0, (4 
=1, t>0. 


Equating these two expressions, multiplying by sin (mzx/a) and integrating over 
(0, a) with respect to x, (m is any positive integer), 


—El,Cns =- Rn (+) sn. 


If the beams of set A are equally spaced, then hx,=ar, where h=p+1, so that 
equation (5) becomes, writing Z,,, for (m=)* Cns, 


Now, for integral values of m and k which are less than h, we have 


r=1 h h 


h-1 

2 sin—-sin—-=h, if n=k (7) 
=0, if nk. 
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Then, from (6), we obtain 


=2 & R,,sin— sin —— 
h h 


This shows that if the Z,,,=(m7)* Cm, are found in the manner outlined in Ref. 1, 
then the reactions R,, at the joints can be found as the sum of a finite series of 
trigonometrical functions, and hence bending moments may be found by taking 
moments. It may be observed that, for relation (8) to hold, the only conditions 
apart from those of the ordinary theory of uniform beams are that one set of beams 
(B) be simply-supported and the other set (A) equally spaced. 


To obtain an expression for direct calculation of deflections or bending 
moments in beams B at the intersections, the procedure is as follows. 


The right hand side of equation (6) is a periodic function of m which is 
unaltered if m be replaced by v= | m+2jh |, (j is any integer positive or negative), 
except for a change of sign if m+2jh < 0. 


Hence Zvs=Zms Sgn (m+ 2jh) (9) 
where sen(j=+1, t>0 
gn (10) 
=-l, t<0 


Then for x=ra/h, equation (1) becomes 


We= & C,, sin—. 
k=1 h 


Rearranging the terms in this series, 


. mtr 
We= C,, sin —— sgn (m+ 2jh). 
m=1 j= h 
Using equation (9), 
We= mCns(m+2jh)-* sin —. 
m=1 j= h 


Now by the known result 


ao - - 
= (m+ 2jh) 48 3 cot ah cosec 


j=-xX 


(11) 
say, 
Tr 


D 
MT 
we obtain Wee Ky, Zune (12) 
m=1 
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Similarly, for bending moment in the B-beam at the joint (r, s) 
M,,=—El, (kr /a)’Cy, sin (kzr/h) 


v*x*C,, sin (mmr/h) sgn (m+ 2jh) 


= m (m+ 2jh) $in (mr /h) 
m=1 j= 
= ms cosec? (mz /2h) sin (mtr /h) (13) 
(14) 


S (m+2jh)-? = cosec” /2h) 


using the summation 


2.2. SECOND CASE: CLAMPED ENDS (SYMMETRICAL IN B-DIRECTION) 


In this case the appropriate series for deflection of the s'" beam of B, assuming 
symmetry about its centre, is 
Cys [1—cos (2kxx/a)] (15) 
k=1 


and, if no external load acts, the bending moment is given by 
D 2 
M,=} [ —(x-—x,)+ (x - x,) U U (x- at x,) | R,, (16) 


Equating this to E/,d*w/dx*, multiplying by cos (2m7x/a), and integrating with 


respect to x over (0, a), 


El, Cu (2) = (1-cos =), 


or for equally spaced A-beams, 


(1-cos 


This shows that Z,,,=0 if m is any multiple of 4, and also that Z,,, remains the 


D 
El, Rr 


same if m be replaced by (h-m). 


Moreover, for 0< m<h and 0O<n<h, 


(1-¢ os =—h if m=n or m+n=h 
j if m=n=4h (18) 
=0 for other values of m and n 
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Hence, from (17), 


EI, = cos 2mrr = S R,. (1-co gmt) eos 


m= h m= 


h 7 


using the property of symmetry, which implies that R,, is unaltered when r is 
replaced by h-r. 


The recurrence relation corresponding to equation (9) is 


where u=|m+jh |, and j is any integer, positive or negative. Thus, proceeding as 
before, the series for deflections at the intersections reduces to 


2mtr 
Wrs >, l — COS 
m 


(20) 


where L,,= (1 +3 cot? =) cos sec* 


and the summation is from 1 to 4p, if p is even, or from 1 to 4(p+1), with the 
final term halved, if p is odd. In either case the expression could be written 


4 
Wrs = — COS - +) 


m=1 h 


Similarly, for bending moments in the B-beams, 


2, mr 
M.,= = Zms COS cosec?—— . 
h h 


3. Solutions for a Particular Grillage 


As an illustration of the method, consider a grillage of three equal uniform 
girders (B-beams) evenly spaced and of length a, supported by nine equal uniform 
stiffeners (A-beams) of length b=4a and moment of inertia one quarter of that of 
the girders (i.e. 4/,=/,). Then p=9 and g=3. We assume a central concentrated 
load P, and consider first all beams simply-supported, and, secondly, all beams 
clamped at the ends. 


With simple-support, we have for the A-beams, by beam theory, 
7b°R,, =6 x (— low,, + 11w,,) 
7b* (R,.— pA,)=6 x VET, (22w,, — 16w,,), 


where A, = 1 and A,=0 for r45. bi 
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Substituting into equation (6), and cer w,, and w,, by their series (12), 


= SK, sin _ 16Z,,, + 11Z,,,) 


=6 x VEI,K,, (- 16Z,,,+ x 10, 


mk 


7b* El 6 x 4°El, sin—~ K, sin 16Z,. 2P sin 


“0 n=! 


=6 x (22Zny — 16Z ns) x 10+ 146*P sin 


Solving these two equations, putting C= 60 x 4° (/,a*)/(1,b*) = 60 x 4 x 3*, 


El, 22PK mC sin (mt /2) 

_2(1+16K,,C) P sin (m=/2) 
a 7432K,C + 2Kn2C? 


After evaluating these expressions for the appropriate values of m, one can 
determine R,, and w,, from (8) and (12), and M,, from (3) for A-beams or from (13) 
for B-beams. Table I gives numerical values for deflection, bending moments and 
shearing forces at the centre of the grillage. The corresponding values obtained 
from a solution by relaxation (taken from an unpublished Admiralty report) are 
also given, for comparison. 


TABLE ! 
NUMERICAL SOLUTIONS FOR THE PARTICULAR GRILLAGE. SIMPLY-SUPPORTED ENDS 

Method of solution Finite series | Relaxation 
Central deflection 7-300 x 10-4 Pa® /(E1,) 7-308 x 10-4 Pa® /(El,) 
Central bending moment: 

(i) Girder (B-beam) — 33-22 x 10-°Pa — 33-02 x 10-4Pa 

(ii) Stiffener (A-beam) —24:12 x 10-5Pa 24°18 x 10-*Pa 
Central shearing force: 

(i) Girder — 28:46 x 10-2P 28°54 x 10-?P 

(ii) Stiffener 10-2P -21-62 x 10-2P 


The analysis for the same grillage with all beams clamped at the ends is as 
follows : — 


Applying beam theory to the A-beams, 
b*R,, =3 x 4°El, (- 8w,, + 4,2) 


b® (R,,—A,P)=3 x 4°El, (8w,, 
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When substituted into (17), w,, being replaced by means of (20), these two 
equations give 


9 
= 3 x [ s > L, (- 8Z,,, + 4Z),.)+ (- mn, +4Zn3)| 
- n=l 


8El,a~*b*Z,,. — (1 —cos mz) 
9 


= 3 x 4E!, | 5 Ly (8Z,., - - 5Zn2) + SLn (8Zn, 


n= 


which may be written 


(1 +8Nn) Zmy — — 8S, + 4S, 
8N Zn, +(1+5Nw) (l— 85, —55, 
= +( mz ~ ( COS Tt) + ) 
where Nn= 120 a®b-* (1,/1,) Ln =810 Ln. 
SEI, (1—cos mz), 
9 
and S.= NiZus 
n=1 
From equations (22), 
8 (14+ N,,) S, +48, (23) 


where An=14+13Nn+8Nm? 


Multiplying these equations by N,,4,,~' and summing with respect to m from 
1 to 9 gives 


[4N - 8 (1+ S,+4S8,] Nm 


9 
> {a + 8N,,) P,.+ 8S, —(5+8L,) N 


m=1 
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The solution of these equations is given by 


AS,=4A,A,+4(14+13A,+8A,) A, (24) 
AS. =(1+8A,+8A.) A,+8(1+12A,+8A.,) A, 
9 9 
where A= 2 Ne Aa = Na? 
m=1 
9 


9 


m=1 m= 


It is then necessary to substitute from (24) into (23) and thence into equation (19), 
(20), or (21) to obtain respectively forces, deflections or bending moments in B-beams 
at the intersections. Numerical results for this case are included in Table II. 


TABLE II 
NUMERICAL SOLUTIONS. CLAMPED ENDS 

Method of solution Finite series Relaxation 
Central deflection 2°674 x 10-4 Pa* /(E/,) 2°678 x 10-4 Pa*® /(EI,) 
Central bending moment: 

(i) Girder — 25:54 x 10-°Pa — 25-52 x 10-3Pa 

(ii) Stiffener - 16°84 x 10-°Pa — 16°87 x 10-8 Pa 
Central shearing force: 

(i) Girder — 26:91 x 10-?2P — 26°891 x 10-?P 

(ii) Stiffener — 23-09 x 10-?P — 23-127 x 10-?P 


4. Conclusion 


As originally presented by Klitchieff, the solution was obtained in the form of 
infinite Fourier series and, for bending moments, direct numerical evaluation was 
tedious if not impossible, owing to the slow convergence of the series. It may be 
seen that computation is made shorter by the modification in the present paper; for 
instance, the example requires the determining of the coefficients Z,,, as far as the 
first five only. 


ACKNOWLEDGMENT 


This paper is published with the approval of the Lords Commissioners of the 
Admiralty, but the responsibility for any statements of fact or opinions expressed 
rests solely with the author. 


REFERENCES 


1. KuitcueFF, J. M. Beams on Elastic Supports and on Cross Girders. The Aeronautical 
Quarterly, Vol. Il, pp. 157-66, November 1950. 


2. KuitcHiEFF, J. M. Beams on Cross Girders with Clamped Ends. The Aeronautical 
Quarterly, Vol. Ill, pp. 230-7, November 1951. 


February 1957 57 


0 
) 


An Approximate Method for the Calculation 
of Turbulent Boundary Layers in Diffusers 


J. F. NORBURY, M.Eng., A.F.R.Ae.S. 


(Department of Fluid Mechanics, University of Liverpool) 


SuMMaARY: An approximate method is described for the calculation of turbulent 
boundary layers in which the turbulence is developed before the commence- 
ment of the adverse pressure gradient, as in most diffuser layers. It is based on 
a method due to Spence®? which has been modified and also extended to the 
calculation of three-dimensional diverging layers such as occur in ducts whose 
breadth is increasing. The velocity profiles occurring in a diverging layer are 
examined and it is shown that the inner part obeys the universal logarithmic 
law, as in two-dimensional layers. This result is used to obtain an equation 
for the form parameter y=(u/U),~» in diverging layers, by substitution in the 
equation of motion and incorporation of the equations of momentum and 
continuity for diverging flow. The form parameter equation contains a term 
involving the gradient of shear stress at y=6 and values of this term are 
obtained by the analysis of experimental data and the substitution of known 
values for all the other terms in the form parameter equation. Values of the 
term involving shear stress gradient are then correlated in terms of known 
boundary layer quantities, and the resulting correlation allows the formulation 
of a step-by-step method for the solution of the form parameter equation. 
This may be used in conjunction with the momentum equation to predict the 
boundary layer growth. It was not found possible to effect a satisfactory 
correlation for boundary layers on lifting aerofoils, in which the turbulence 
develops within the adverse pressure gradient, and the method cannot be used 
for the prediction of such layers. The results of a number of calculations 
are given. 


1. Introduction 


In recent years a number of investigators have considered the problem of 
calculating the growth of turbulent boundary layers. The mechanism of such a 
layer is very complex and is not by any means fully understood. Thus a completely 
theoretical treatment is not at present possible, and the semi-empirical approach 
is difficult. 


The problem is complicated by the variety of conditions which may occur at 
the beginning of the turbulent layer. In boundary layers on aerofoils, the transi- 
tion to turbulence may follow an instability of the laminar layer or a bubble of 
separated flow. The development of the turbulence is presumably dependent on 
the severity of the adverse pressure gradient in the transition region, and it seems 
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possible that the nature of the transition influences the subsequent growth of the 
layer for a considerable distance. In diffuser flows the boundary layer is usually 
turbulent at the throat, but its thickness and velocity profile, and the state of the 
turbulent motion, depend on upstream conditions and may vary considerably. 


Usually the object of a boundary layer calculation is the determination of the 
momentum thickness @ and the velocity profile form parameter H, the product H@ 
giving the displacement thickness 5*. The principal difficulty lies in the determina- 
tion of H. For typical applications this is not required to a high degree of accuracy, 
and any method which would unfailingly predict H within five per cent would 
probably be adequate. Such an accuracy cannot yet be attained. 


The present work arose from a desire to analyse boundary layers occurring in 
a diffuser of rectangular section and to assess the possibility of calculating the 
growth of the boundary layers along the centre line of each wall. The roof and 
floor of the diffuser were plane and parallel, and the side walls were flexible and 
could be moved to vary the pressure gradient. Thus the boundary layers on the 
side walls might be expected to be two-dimensional over the early part of the 
diffuser, where all the layers were relatively thin, and the flow in the boundary 
layers on the roof and the floor was divergent. 


Two recent methods for the calculation of turbulent boundary layers were 
considered. That due to Rubert and Persh") depends on the kinetic energy 
equation and the assumption of power law velocity profiles for the formulation of 
an equation for H. It was applied in the original paper to the calculation of both 
two-dimensional layers and three-dimensional layers occurring in conical diffusers. 
However, calculations by this method did not correspond very well with the 
experimental measurements in the present case, the disparity being especially 
marked in the case of the three-dimensional layers. The method due to Spence”? is 
based on the proposition that the velocity profiles van be expressed by means of 
the logarithmic law over the inner part and a power law over the outer part. This 
appears to be fairly satisfactory and the method is convenient to use, but in its 
original form there are features which render it unsuitable for calculations of diffuser 
boundary layers. Modifications are here presented to avoid these difficulties, and 
the method has been extended to the calculation of the three-dimensional diverging 
boundary layers occurring in the diffuser. 


NOTATION 

x,y co-ordinates parallel and normal to bounding plane of flow 

u,v mean velocity components in x, y, -directions 

U _ free stream velocity at edge of boundary layer 
p,¥ Viscosity, density, kinematic viscosity 

tT, Shear stress, shear stress at wall 
c: local skin friction coefficient =-,,/(4pU’) 
u; friction velocity=(7w/p)!=U (c;/2)! 
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8 thickness of boundary layer 


5* -- displacement thickness of boundary layer 


II 


“(1 “Yay, momentum thickness of boundary layer 


6* = 

H = 8*/6 
V>= (u/U)y=6 

R, = U6@/», Reynolds number based on momentum thickness 


y* = yur/¥ 
a,,k, constants occurring in Spence’s method 


} orm parameters for velocity profile 


x(y)_ function occurring in Spence’s method 


X distance expressing the divergence rate of a three-dimensional 
diverging layer 


f function for inner part of velocity profile 


F,2,G,,G, functions occurring in form parameter equation and defined in 
equations (12), (14), and (16) 


k quantity defined in equation (11) 
U =U/U, 


Suffixes 
refers to initial values 


refers to three-dimensional diverging layer 


2. Outline of Spence’s Method 


It is required to calculate the variation of 0, H, 5* and 7 along the length of 
the boundary layer. The four equations used are 


(i) the boundary layer momentum equation, 
(ii) the relationship H =6*/6, 
(iii) an empirical equation for /(4pU*), 


(iv) an equation for the variation of the form parameter H derived from the 
equation of motion. 


To obtain the fourth equation, Spence’ introduced as basic form parameter 
the quantity y defined as the ratio of velocity at distance @ from the wall to the 
main stream velocity U. 
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For a power law profile 


H-1 . 
and, as this relation is nearly correct for experimentally measured layers, it is used 
in the subsequent work. A differential equation for the variation of y along the 
layer was obtained by substitution in the equation of motion of the logarithmic law 
for velocity distribution near the wall, 


[2:5 log.( Uy +5:5| @ 


This differential equation involved the gradient of the shear stress at a distance 
from the wall, (07/dy),-2. No theory exists by which this quantity might be 
predicted and it is difficult to determine experimentally. Spence introduced the 
relation 


and then determined the constants z, and k, so that his form parameter equation 
gave a good fit with experimental results. He then reduced the form parameter 

n equation to a more convenient form in which it may be used in conjunction with 
the momentum equation to calculate the boundary layer growth by an 
iterative process. 


In attempting to apply the method to the calculation of boundary layers in 
diffusers, difficulties arise on two points, namely the assumptions regarding the 
shear stress gradient and the choice of the initial value of H. In the throat of a 
diffuser U’/U =0 and, in general, 


dy pU? v=0 


Thus 2, becomes infinite and, in the neighbourhood of the throat, where U’/U is 
small, <, will still be considerably different from any constant value which may be 
suitable over the greater part of the layer. The choice of a suitable initial value 
for H is important for satisfactory calculations. For aerofoils, Spence indicated 
that a constant initial value, actually 1:4, gave good results in the calculation of 
layers having free transition. For a diffuser, the initial value of H may vary between 
wide limits, and it is not possible to assume a fixed value. If the initial value of H 
is very low the Spence method is inapplicable. The function 


x ()=1-797 (122+ 


r which appears in the final form of Spence’s form parameter equation actually has 
; no real values for y > 0-82, corresponding to H < 1:14. In fact values of H less 
than 1-14 may occur at diffuser throats and such cases could not be dealt with. 
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Fig. 1. 
Velocity profiles in diverging flow plotted logarithmically. 


3. Velocity Profiles in Diverging Boundary Layers 


For two-dimensional layers, Spence showed that in an adverse gradient the 
inner part of the velocity profiles could be expressed by the logarithmic law, 
equation (2), and further that the values of the parameter (u,/U) required for 
correspondence were approximately equal to the values of (c;/2)! (=u,/U by 
definition) given by the skin friction relations of Ludwieg and Tillmann™? and, for 
low values of H, of Squire and Young™’. These relations are, respectively, 


and c;=2 (3:59 + 2:56 log. R»)~* ‘ ‘ (6) 


A similar investigation has been made for a three-dimensional diverging layer 
occurring in a diffuser, the velocity profiles having been measured along the axis of 
symmetry of the layer. The validity of the logarithmic law over the inner part. 
approximately one-fifth of the layer, is shown by the results plotted in Fig. 1. No 
reliable measurements of the skin friction coefficients for the layer were made, so 
that it was not possible to make a direct check of the accuracy of equation (2). The 
course adopted was to assume that equation (2) was valid for the three-dimensional 
flow, and then to compare the resulting values of (u-/U) with (c;/2)! given by the 
skin friction relations for two-dimensional flows with equal values of R, and H. 
In fact the results showed a close resemblance to those obtained for the two- 
dimensional layers. The values of u-/U indicated by the velocity profiles fell 
slightly below those given by the Squire-Young relation (Fig. 2). For the 
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Fig. 2. 
Variation of skin friction in diverging flow. 


Ludwieg-Tillmann relation the disparity was greater over the early part of the 
layer, where the values of H were rather low, being about 1-25. 


In addition the relation between y and H for the diverging layers corresponded 
closely with that given by the power law, the order of disparity being the same as 
that for two-dimensional layers. 


4. Form Parameter Equation for Three-Dimensional Diverging Flow 


We are considering flow over a plane boundary in which the streamlines 
diverge away from some central axis of symmetry, and we require an equation 
expressing the variation of the form parameter equation along the axis of symmetry. 
This is taken to be the x-axis, the y-axis being at right angles to the plane boundary. 
u and wv are the corresponding velocity components. The streamline divergence 
rate is conveniently expressed by the distance, X, from the section considered to the 
point of intersection of the x-axis and a tangent to the streamline at the section 
(Fig. 3). In general the streamlines are curved, so that X will vary along the length 
of the layer, and at any particular section it may also have different values for 
streamlines at different distances from the axis of symmetry. The required value 
of X is the value which is approached for streamlines approaching the axis 
of symmetry. 


In the laminar sub-layer, 
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Fig. 3. 
Determination of X from streamline pattern. 


and just outside this is a region in which the logarithmic law, equation (2), is valid. 
The momentum thickness @ lies within this region (Fig. 1). 


Thus near the wall, 


Vu; 
— =f (y*) where y*=-—., 
u; 


The continuity equation for the type of flow considered is 


du ov eu 


while the equation of motion in the direction of the x-axis for steady flow is the 


same as for two-dimensional flow, 


Ou aU ldr 


The momentum equation becomes 


Writing out the equation of motion at y=, and using equations (2), (7), and (9), 
we obtain 


ur Ur u;\? F (6*) 
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Regarding y*=11-7 as the outer limit of the laminar sub-layer, we find 


F (0) (8) -6-25- (13) 
ur Rott 
Introducing also the function 


equation (10) becomes 


where G,= Sy EY) and (16) 


For two-dimensional layers X = 090 and the form parameter equation becomes 


Curves of the functions g, G, and G, for various values of H and k are given in 
Figs, 6, 7, and 8. Tables of values have been calculated but are not presented here. 


5. Calculation of Two-Dimensional Layers 
5.1. INTRODUCTION AND CORRELATION OF EXPERIMENTAL DATA 


Since the product (u,/U)(H+2) changes little over most of the layer, the 
variation of the quantity k, defined in equation (11), depends largely on the shear 
stress gradient at y=. The variation of k for a number of experimental layers 
was found by the substitution of experimental data in the two-dimensional form of 
equation (10), the Squire-Young relation being used for the calculation of u-/U. 
The following series of layers was analysed : — 


(a) Schubauer and Klebanoff’’, 


(6) McCullough and Gault. Boundary layers on N.A.C.A. 63,-018 at 0°, 
6°, 12° incidence, and on N.A.C.A. 63,-012 at 12-8°, 


1) (c) Altman and Hayter’. Unswept wing at C,=0, 0-32 and 1-0. 


(d) von Doenhoff and Tetervin’’. Nose-opening aerofoil, shape 13, 91° 
incidence, and N.A.C.A. 65 (216)-222 at 10°1°. 


(e) Measurements made by the author in a diffuser of rectangular section. 
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Fig 4. 
Values of k from experimental data. 


An attempt was made to correlate the values of k which were obtained, in terms 
of known boundary layer quantities. After preliminary trials the two parameters 
chosen for the correlation were R, and H,* 6U’/U. 


(H, is the initial value of H. For those layers in which transition was complete 
before the adverse pressure gradient began, the initial value was taken to be that at 
the minimum pressure point. For the layers on aerofoils at incidence, where the 
end of the transition region was not clearly defined, the initial value was taken to 
be that value for which dH/dx=0). Values of k were plotted against R», and the 
resulting distribution of points is shown in Fig. 4, where also the value of 
H,° (@U’/U) for each point is given. Negative values of k were obtained in some 
cases from points near the start of the adverse pressure gradient. It may be 
supposed that at such sections the maximum shear stress lies closer to the wall than 
the distance y=9, so that (07/d0y),-s is negative. 


It appeared that the manner in which k varied was different for two distinct 
groups of layers; namely, those in which transition was probably complete before 
the commencement of the adverse pressure gradient (the Schubauer-Klebanoff layer. 
the diffuser layers, and those on aerofoils at zero incidence), and those in which 
transition occurred within the adverse pressure gradient, as on the aerofoils at 
incidence. Points obtained from the latter group could not be correlated by the 
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Fig. 5. 
k—correlation chart. 


parameters used and this applied to points derived from the downstream parts of the 
layers as well as to those from the parts nearer to the minimum pressure position. 
It is possible that the lack of correlation of the values of k from such layers reflects 
the complex nature of events in the transition region, and it suggests that the nature 
of the transition influences the whole subsequent behaviour of the turbulent layer. 


For the diffuser layers, on the other hand, the values of k were correlated 
reasonably well. It was possible to draw lines of constant H,* (6U’/U) on the 
diagram of k against R», and the resulting chart (Fig. 5) could then be used to 
determine k for boundary layer calculations. It is difficult to understand the 
physical significance of the parameter H,* (6U’/U), and the correlation of k values 
must be regarded as purely empirical. 


5.2. METHOD OF CALCULATION 


Before calculating the variation of the form parameter the momentum 
equation must first be solved, using H=constant, H,, to obtain a first approxima- 
tion for 6. This may conveniently be carried out by the method described by 
Spence’), The quantities 4, U?, U’/U, Re, (1/9) (uz/U) and 6U’/U are then 
obtained for all points of the layer. The terms Ro, H,,* (@U’/U), (1/9) (u-/U)*, and 
U’/U may then be plotted and equation (17) may be solved step by step. The term 
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Fig. 6. 
Curves of function g (k, H). 


on the right of the equation is small by comparison with U2g, and over the greater 
part of the boundary layer large steps may be taken. The size of the step allowable 
may be determined by the rate at which k varies, if this is rapid. Over most of the 
layer k varies only slowly. 


The procedure to be followed for each step is: — 


(i) Choose size of step, Ax. Determine mid-step values of R, and 
(6U’/U). 


(ii) Select k from chart (Fig. 5). 

(iii) Knowing k and H, find g from curves (Fig. 6) and hence Ug at beginning 
of step. 

(iv) Find mid-step values of (1/6)(u-/U)’ and U’/U. Find G, and G, from 
curves (Figs. 7 and 8). Hence find o (U?g) for the step. 

(v) = (U?g) Ax. 

(vi) Find U?g and hence g and H at end of step. 


When the form parameter equation has been solved, a second solution of the 
momentum equation, using the known variation of H, may be obtained, again using 
Spence’s method’. This will usually complete the boundary layer calculation to a 
sufficient degree of accuracy, but further solutions for H and @ may be obtained 
if necessary. 
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Fig. 8 
Curves of function G, (k, H). 
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Throughout the calculation it is convenient to use the Squire and Young 
expression for skin friction, rather than that of Ludwieg and Tillmann. For the 
solution of the momentum equation the difference between the two relations is small, 
except near separation. when the skin friction term is no longer very important. 
In any case the higher value given by Squire and Young helps to offset the omission, 
in the usual form of the momentum equation, of terms arising from the boundary 
layer turbulence. 


Two further points arising in the derivation and use of the method remain to 


be dealt with. In obtaining the derivative £ (Ug) it was assumed that k was 


independent of x. This condition is satisfied for each individual step. 
The function 2 (k, . (14 


becomes zero when y= 1—k, and for y > 1-k it has no real values. In practice it 
was found that the use of values of g in the neighbourhood of zero gave unreliable 
results, and on the g, H-chart (Fig. 6) a line was drawn to indicate the boundary of 
this region. If the k-correlation chart indicates the use of a (g, H) point falling 
below the line, the procedure is to choose the next lower value of k for which a curve 
is available, and for which the initial value of g for the step falls above the line. This 
procedure is not often necessary and must be regarded as a secondary correlation 
for k which over-rides the normal one for low values of H. 


5.3. DISCUSSION OF CALCULATIONS AND RESULTS 


The results of some typical calculations of H are shown in Figs. 9-14. The 
several quantities required in the calculation were taken directly from experimental 
data and not from a first solution of the momentum equation. Not more than ten 
steps were used for any of the calculations. 


The layer on a symmetrical aerofoil at zero incidence was calculated by the 
methods of Spence and of Rubert and Persh, as well as by the present method, and 
the results are compared in Fig.9. There is considerable scatter of the experimental 
points and it is not easy to choose a starting point. Curves from H,= 1:34 and 1-35 
calculated by the present method enclose a mean path through the points, and 
indicate the influence of a small variation in H,. The Rubert and Persh result lies 
appreciably below the experimental points although it was begun at a higher initial 
value. The negative slope of the Spence result over the early part of the layer is 
indicative of the difficulty occurring in this method in a region of low 
pressure gradient. 


Figure 10 shows a calculation for a layer on one wall of a rectangular section 
diffuser. Over the early part, where the boundary layer is thin and may be 
expected to be two-dimensional, the correspondence is good. At about half the 
distance along the layer the effect of interference of the boundary layers on the 
adjacent walls is apparent and this influence increases steadily. In Fig. 11 is shown 
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Growth of H on N.A.C.A. 63,-018 at 0° (Ref. 6). 
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Fig. 10. 
Growth of H on wall of rectangular section diffuser. 
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Fig. 11. 
Growth of H on wall of rectangular section diffuser. Layer artificially thickened 
upstream of throat. 


the result of a calculation of a similar layer which was artificially thickened 
upstream of the diffuser throat so that it had a very low initial value of H. In 
calculating this layer the k value had to be modified for several steps to avoid the 
invalid region of the g, H-chart (Fig. 6). However, the calculation is successful up 
to the point at which three-dimensional effects become apparent. 


The result for the Schubauer-Klebanoff layer is shown in Fig. 12. Since this 
experimental layer completely determined the form of the k-correlation chart for 
Rs > 20,000, the calculation must be expected to yield a good result. Calculations 
were made from two starting points, H,=1-:30 and 1-35, to illustrate the influence 
of the initial H value on the result for a separating layer. Clearly an initial point 
can be chosen somewhere between these limits to give very good correspondence 
between calculated and experimental results. 


The result of a calculation for a layer on a plane wall investigated by 
Sandborn™? is shown in Fig. 13 and compared with a calculation by the method of 
Maskell". It is interesting that the two calculations correspond so closely, 
although the correspondence with the experimental results is not so close. Actually 
the boundary layer occupied a large part of the depth of the duct, and there is 
reason to doubt whether it was truly two-dimensional. 
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Growth of H for Sandborn layer (25 in. suction). 
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6. Calculation of Three-Dimensional Diverging Boundary Layers 
6.1. CORRELATION OF EXPERIMENTAL DATA 


For the purpose of analysing the diverging boundary layers occurring in the 
rectangular diffuser, the divergence rate, expressed by the distance X, was found by 
plotting the streamline pattern for potential flow using an electrical method. Values 
of ka for the layers were then worked out from equation (10). (The suffix “d” 
distinguishes values of k relating to diverging layers). These values were slightly 
but systematically higher than the values indicated by the k-chart for two- 
dimensional layers having the same values of R, and H,* (6U’/U). The differences 
(ka-k), were plotted against corresponding values of 6/X (Fig. 14). A roughly 
linear relationship was indicated, points being scattered about the line 


Since the maximum difference was only 0:1, this line was quite adequate to 
determine the correction to be applied to the two-dimensional value for the purpose 
of calculating three-dimensional layers. 


6.2. METHOD OF CALCULATION AND RESULTS 


The form parameter equation to be solved is 


4Up=0'| -G, + {625+ (15) 
U 
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Fig. 15. 
Growth of H for diverging layer in rectangular section diffuser. 


For all but small values of 6, 96/(Reu-/U) is much smaller than 6°25, and when @ is 
small the third term on the right hand side is much less important than the first, 
which itself is small compared with U*g. Thus the calculation may be simplified 
without loss of accuracy by writing the third term 


G. 1 (u-\? 
= 


G,/k,? can be determined from the curves for G. and the known value of ka. It 
! must be noted, however. that as ky approaches zero, G, approaches zero, while 
G,/k,? approaches 2. 


Apart from taking account of this third term, the calculation differs from that 
for two-dimensional layers, described in Section 5.2, in only one feature. In step 
(ii), after selecting k from the k-chart (Fig. 5), the correction for divergence, 
+55 6/X, must be applied. Since the primary correlation of k is slightly different 
for the two-dimensional and divergent boundary layers, it may be expected that the 
secondary correlation, represented by the limitation on the valid region of the 
ge, H-chart (Fig. 6), will also differ slightly. However the use of the two-dimensional 
limit line gives satisfactory results, provided that the variation in k which it 
dictates is less than 0-1. 


The result of a calculation, using experimental data directly, for one of the 
diffuser boundary layers is shown in Fig. 15, together with that of a calculation for 
the same layer by the method of Rubert and Persh. 
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7. Conclusions 


The analysis of data for the determination of values of k indicated a consider- 
able difference between turbulent boundary layers on lifting aerofoils and those 
occurring in diffusers. No correlation of data for boundary layers of the former 
type could be obtained by the method used. This may indicate the influence of the 
complex mechanism of transition, and may suggest that the uncertainty surround- 
ing the transition region of aerofoil layers is to some extent responsible for the 
formidable difficulties encountered in attempting to predict the growth of 
such layers. 


For the diffuser boundary layers, a correlation of data was possible and from 
this a method of calculating the growth of the form parameter H has been evolved. 
This has given reasonable results for experimental layers from several different 
sources. Although 2 step-by-step calculation is required, this is not particularly 
laborious and it appears that as much numerical accuracy as the method justifies 
can be attained by the use of a relatively small number of steps, perhaps about ten. 
The application of the methods to diverging layers is only slightly more complicated 
than the application to two-dimensional layers. 


The formulation of the method required the introduction of experimental 
evidence concerning the distribution of shear stress across the layer, this being 
embodied in the correlation of values of k. This correlation was purely empirical 
and cannot be regarded as completely satisfactory. An important requirement in 
the field of turbulent boundary layer calculation is a theory for the prediction of 
the distribution of turbulent shear stress. 
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Stability and Response of Systems Satisfying 
a Second-Order Linear Differential Equation 
with Time-Dependent Coefficients 


A. W. BABISTER 


(The University, Glasgow) 
SUMMARY: The differential equation considered is 
d*x dx 
(a, + +(a,+,7) +(a,+b,7) x=0, 


where all the a’s and 5’s are real constants. 


The nature of the solution is investigated in the neighbourhood of the 
singular point (r= — a,/b,), and the conditions are found for logarithmic terms 
to be absent. 


The conditions for stability for large values of 7 are determined; the system 
is stable if 


{a,+5,7), (a,+6,7) and (a,+5,7) 
are all positive for large values of 7. 


The form of the response is considered and its oscillatory (or non-oscilla- 
tory) naiure investigated. The Sonin-Pdélya theorem is used to determine simple 
inequalities which must hold between the coefficients of the differential equation 
in any interval for the relative maxima of |x| to form an increasing or 

decreasing sequence in that interval. 


I. Introduction 


A system is considered which satisfies the linear equation 


(a, + b.7) ax +(a,+b,7)+(a,+6,7)x=0 . (1) 
dr? dz 


where all the a’s and b’s are real constants and b, ~ 0. An equation of this type can 
arise (i) where a linear system is operated over a sufficiently wide range, (ii) for 
mechanical systems in which mass, damping or stiffness varies with time and (iii) 
for electrical systems operating while warming up in which inertia, damping or gain 
varies with time (see also Ref. 1). 
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The initial conditions are given by 


The general solution of (1) is 


where C, and C, are constants satisfying the initial conditions (2) and x, and x, are 
any two independent solutions of (1). In general x, and x, can be expressed in 
terms of confluent hypergeometric functions (see Refs. 1 and 2). Some general 
properties of solutions of (1) are now considered. 


In general, b, can be taken to be positive without loss of generality. We are 
only interested in values of x for positive values of 7. 


2. Singular Points 


Equation (1) has one regular singularity at >= - a,/b,. We write 


If b, is zero there is no singularity at a finite time. If a, and b, have the same 
sign (positive), the singularity will occur for a negative value of = and is thus not of 
physical importance. If @, and b, have opposite signs it is found that, in the 
neighbourhood of this singular point, the solution of (1) can be put in the form 


ZA, x, (r), 
33 


r 


where x, (7) is given by 


|. 


The indicial equation for 2, is 


2 a,b, a,b, 


The roots of the indicial equation (5) are 0 and 1-c. Thus the solution of the 
differential equation corresponding to the zero root of the indicial equation is always 
finite at the singular point. 


If c > 1, one solution of the differential equation will tend to infinity at z,. If 
¢ < 1, both solutions of the differential equation are finite at ,, one solution tending 
to zero as 7 approaches ;,. If 7, >0, this solution is in general indeterminate 
for r > 1,. 
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If c is an integer, one solution of the differential equation in general contains 
terms such as 


(7 log (7 = To). 
and there is a logarithmic singularity at ;,. For c< 1, the exceptional cases in 


which the logarithmic terms are absent (and both solutions are finite at z,) can be 
found as in Ref. 3. 


ay 
For simplicity we put 


Then (1) becomes 


ax dx 


Sig +(A,+ Bis) 1, +(A,+B,s)x=0 (7) 
where A,= = 
b, 
= b, 
b (8) 
A, b.? 
by 
B,= b, 
co 
Let x= ‘ : (9) 
v=0 
Substituting in (7) and equating coefficients, we have 
Coefficient of 
+2) (6 + 1)+ A, +2)] [B, (6 +1)+A,]+c,B,=0 . (12) 


ses +3) (6 +2)+A, (6 +3)] +, [B, (6 . (13) 


As shown, the roots of the indicial equation (10) are 0 and 1-A,. Consider 
negative integer values of A,; the two solutions of the differential equation will 
contain no logarithmic terms if C, a, iS indeterminate when we put o=0 in 
equations (10)-(13). 
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Putting «=0, equations (11)-(13) become 
c,A,+c,A,=0 
c,(2+2A,)+c, (B, +A,)+c,B,=0 
(6+3A,)+c, (2B,+A,)+c,B,=0. 


If A, <1, from these equations it is found that the following relations must hold 
between the coefficients for logarithmic terms to be absent. 


c=A,=0, A,=0 
c=A,=-1, |A, -1 

|B,  B,+A,| 

ta 0| =0 

| 

| B, B, + Ay | 

0 B, 2B,+A, | 


3. Singular Point at the Origin 


When a, is zero, the singular point occurs atr=0. From (9) and (10), it is seen 
that one solution of the differential equation will make x, infinite if c>0 
and x, infinite if c> -1. Similarly a solution containing logarithmic terms will 
make x, and x, infinite. 


If c< -1 and is not an integer, both solutions of the differential equation 
make ¥, finite. Then, from (1) with a,=0, 


In this case we are not free to choose both x, and x,. Two independent 
solutions can be obtained by specifying, for example, the values of x at the origin 
and at one subsequent time. Equation (15) also holds when c is a negative integer 
if the particular relations given by (14) are satisfied. 


As already stated, whatever the value of c, there is always one solution of the 
differential equation which is finite with finite derivatives at the singular point. 
However, when we are restricted to this one solution, x, and x, are related by (15). 


4. Stability 


The system satisfying equation (1) is said to be stable if, following any initial 
arbitrary disturbance in position x, or in velocity x,, the system ultimately returns to 
its equilibrium position in which both x and i are zero. If the system only returns 
to its original position for particular sets of values of (x,, x,) the stability is said 
to be conditional. 


Consider the motion at any finite time; x may diverge to infinity only at the 
Singular point r=r,. As already shown, this occurs if c=1. We are concerned 
only with those motions which persist at times subsequent to -,. 
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As shown in Ref. 3, the solution of (1) for large = is given by 
(a, + + ‘ J . (16) 


where 4, is a root of the equation 


B(A)=b,r7 +b,A+b,=0. ‘ . (17) 

A (A) 


If all the roots of (17) are negative or have negative real parts, i.e. if with 
b, > 0, we also have b, > 0 and b, > 0, the system is stable for large values of -. 
It is seen that, as would be expected, if neither b, nor b. vanishes, the conditions for 
stability for large > are identical with those for the motion defined by the equation 


d°x dx 


b, +b, 


If b, is zero, one root of (17) is zero and the corresponding value of ¢, is 


For stability for large 7, with b, > 0 and b,=0, we must have b, > 0 and a, > 0. 


If b, is zero and b, and b, have the same sign, the roots of (17) are of the 
form +iw, where 


The corresponding values of c, are 


_ 4, i 
2b, ~ 2b, 
For stability for large +, with b, > 0 and b,=0, we must have b, > 0 and a, > 0. 


If b, is zero, equation (1) is 


d°x dx 
+ (a, + b,7) +t (a, +by7) x=0. ‘ 


The stability of this system can be deduced from its exact solution. 
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the differential equation for u is of the form 


@u 


+(A + Br+Cr*)u=0. 


This can be solved in terms of the parabolic cylinder functions». It is found that 
for stability with large 7, with a, > 0 and 5, and b, not zero, we must have b, > 0 
and b, > 0. 


If both b, and b, are zero, but b, is not zero, we find that for stability with large 
7 with a, > 0 we must have a, > 0 and b, > 0. 


If both b, and b, are zero but b, is not zero, we find that for stability with large 
7 with a, > 0 we must have b, > 0 and a, > 0. 


Finally if b,=b,=b,=0, the system is stable if a, > 0, a, >0 and a, > 0. 
All the above conditions for stability for large values of = can be expressed by 
the inequalities 


(a, +b,7) > 0, (a,+6,7)>0 and (a,+6.7) >0 


for large values of 7. 


5. Oscillatory or Non-Oscillatory Response 


Using the Sturmian theory” we shall derive some general properties of the 
response of a system satisfying the equation (1) following some initial disturbance in 
position or velocity. 


As before, for simplicity we put 


= + a, =T=- 
To 
Then (1) becomes 
d*x dx 


where A,, B,, A, and B, are given by (8). 
Equation (23) can be put in the form 


d ( 8s A,dx Bs A, 
e 


ds ds Ss (A, + B,s) x=0 (24) 


Only positive values of s are considered. 
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Using Rolle’s theorem” it is seen from equation (24) that between two 
consecutive positive zeros of dx/ds there lies at least one real zero of (A, +B,s) x. 
But between consecutive real zeros of x, by Rolle’s theorem, there lies at least one 
real zero of dx/ds. 


Thus, in any positive interval in s which does not include s=- A,/B, (ie. 
7=-a,/b,), between successive maxima and minima of x there is always one 
real zero. 

A 
Let P=e's 
(25) 


B,s A 


(A, +B,s) 


Then for any finite interval (a, b) such that 0 <<a<s <b, P and Q are continuous 
real functions of the real variable s. P is positive and dP/ds is continuous through- 
out this interval. 


Using Sturm’s fundamental theorem, it is found that if the solutions of (23) 
oscillate in the interval (a, b) they will oscillate more rapidly if P is decreased or 
Q is increased. Thus the response is made less oscillatory by decreasing A, or B, 
(for s > 0). 


Let P;, and P;, Qu and Q, be the upper and lower bounds of P and Q in the 
interval (a, b) already defined. Comparing (24) with 


d dx 


it is seen that the solutions of (23) do not oscillate more rapidly than those of (26). 
If Oy is negative, both solutions of (26) are non-oscillatory. Thus if Q is negative, 
i.e. if (A, +B,s) is negative throughout the interval (a, b), the solutions of (23) are 
non-oscillatory in that interval. 


Comparing (24) with 
d 
(Pe *) +Q x= 0 (27) 
the solutions of (23) oscillate at least as rapidly as those of (27). If Q, is negative, 


the solutions of (27) have a period 27/ (Py/Q;). 


If A,, B, and B, are positive and A, > 1, both P and Q increase monotonically 
for all positive s, and thus the upper and lower bounds of P and Q will occur at the 
ends of the interval (a, 5). 


If b,=0, equation (1) becomes 


d*x 


a, de +(a, +b,7) + (a, +b,7) x=0 (28) 
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Equation (28) can be put in the form 


5 & a, 2a ) = 
*) + x=0. (29) 


This can be analysed as before. It is found that the response is made less oscillatory 
by decreasing a,/a, or b,/a,. If (a,+,7) is negative throughout the interval (a, b) 
the solutions of (28) are non-oscillatory in that interval. 


6. Damping of the Response 


The Sonin-Pélya theorem” can be used to determine whether the successive 
maxima and minima of any solution of (1) form an _ increasing or 
decreasing sequence. 


Expressing (1) in the form (24), it is seen from (25) that if (A, + B,s) is positive 
through the interval (a, b) defined in Section 5, both P and Q will be positive and 
have continuous derivatives. From the Sonin-Pélya theorem, the relative maxima 
of x form an increasing or decreasing sequence according to whether PQ is decreas- 
ing or increasing. 


2B,s  2A,- 


Now PQ=e (A, +B,8) 
19. (A, + B,s) (A, + B,S)- Ac] . 


Thus if (A, + B,s) is positive throughout the interval (a, b), the relative maxima 
of | x |will form an increasing or decreasing sequence according to whether 


(A,+B,s)(A,+B,s) < or > 4A, . (30) 


throughout the same interval. 


Certain particular cases are of interest. If A,, A,, B, and B, are all positive, 
the relative maxima of | x | will decrease for all positive s if A, >4. If B, is zero 
and A, is positive, from (30) the relative maxima of | x! will form an increasing or 
decreasing sequence according to whether 


Bs < or > }-A, ‘ GI 


If B, is zero and (A,, + B,s) positive, from (30) the relative maxima of | x | will form 
an increasing or decreasing sequence according to whether 
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If b,=0, equation (1) can be put in the form (29). This can be analysed as 
before. It is found that with a, > 0, if (a,+ 5,7) is positive throughout the interval 
(a, b) the relative maxima of | x| will form an increasing or decreasing sequence 
according to whether 


< or > O. . (3) 
throughout the same interval. Thus if b,=0 and a), a,, a., b, and b, are all positive, 
the relative maxima of | x | will decrease for all positive 7. 


If both b, and b, are zero, from (33) with a, >0 and (a,+ 5,7) positive, the 
relative maxima of | x| will form an increasing or decreasing sequence according 
to whether 


< or > -(a,a,+4a,b,) . (34) 


If both b, and b, are zero, from (33) with a, >0O and a, >0, the relative 
maxima of | x | will form an increasing or decreasing sequence according to whether 


Finally if b,=b,=b,=0, a, > 0 and a, > 0, the relative maxima will increase 
or decrease according to whether 


se < 
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Note on the Extension of Evvard’s Method 
to Wings with Subsonic Leading Edges Moving 
at Supersonic Speeds 


G. J. HANCOCK 


(Queen Mary College, University of London) 


SuMMARY: Evvard’s technique is applied to the problem of a thin finite wing 
moving at a supersonic speed when the leading edge is subsonic. It is developed 
in two methods : — 


(i) in which the relationship between the pressure loading and the 
integrals of the downwash over the wing surface is extended as far 
as possible, and which has to be computed numerically; 


(ii) in which approximations are made for the upwash velocities in the 
neighbourhood of the leading edge, resulting in a series of standard 
integrals for the estimation of the pressure loading. 


Method (ii) is applied to the pressure loading on a flat plate triangular wing 
and cropped delta wing, and the application to more general shapes is discussed. 


1. Introduction 


The general problem of the determination of the pressure loading on a thin 
wing moving at supersonic speeds when the leading edge is supersonic can usually 
be solved, within the framework of the linearised theory, by Evvard’s method"’. In 
this method the pressure at a point on the wing is given as an integral of the down- 
wash over certain regions of the wing surface upstream of the reference point. 
However, the problem for a thin wing with a subsonic leading edge has not been 
investigated fully from this point of view, although particular cases, such as 
triangular wings and cropped delta wings, have been investigated by alternative 
techniques (e.g. harmonic transformation theory, operational methods’, doublet 
line considerations” and the superposition of conical fields’). 


It seems worthwhile therefore to determine the possibility of extending Evvard’s 
method to this problem of thin wing with a subsonic leading edge because : — 


(i) Its generality covers plan form and camber, and is particularly useful in 
aeroelastic problems. 
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(ii). It is generally recognised that the transonic region is a region in which 
oscillatory troubles are apparently concentrated (e.g. flutter, negative aero- 
dynamic damping) so that oscillatory derivatives are needed. At present 
these are mainly derived from slender-body theory, which is inadequate 
for a full investigation. However, since it has finally been shown that 
Evvard’s method can be applied to the problem of an oscillatory 
rectangular wing’ ”), a path is indicated for obtaining similar derivatives 
of subsonic leading edge wings once the steady case has been investigated. 


The first attempt on this problem was made by Evvard") where, as a first 
approximation to a perturbation velocity potential at a point P on the wing, he 
integrates the downwash over the Mach reflection “ diamond ” § formed upstream 
of the reference point P, shown in Fig. 1. This method can only be applied with 
assurance when the wing leading edge is only slightly subsonic. The error, as 
compared to the exact linearised solution, increases as the Mach number approaches 
unity, to the order of 30 per cent at M=1 for a triangular plate aerofoil. A second 
investigation has recently been published by Etkin and Woodward’, where they 
state that a second approximation to the problem is the inclusion of a further area 
of integration with Evvard’s area S, namely S,, indicated in Fig. 1. 


This approximation is shown to be nearer to the accurate linearised pressure 
distribution over a triangular plate than Evvard’s, with a maximum error of 12 per 
cent at M=1; however, no attempt is made apparently to show that this approxima- 
tion is valid for all wings or to indicate the way to the determination of a more 
accurate solution. 


, MAIN SUPERSONIC 
STREAM 


MACH LINES 


LIFTING 
SURFACE 


Fig. 1. 
Regions of integration suggested by Evvard, and by Etkin and Woodward. 
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In this paper, a complete solution is not derived for the general problem of 
establishing a relationship between the pressure loading and integrals of the down- 
wash over regions of the wing surface. Therefore, two courses are open for 
investigation, as in subsonic incompressible finite wing theory. Either, integral 
formulae must be developed which can be dealt with numerically, so that the 
accuracy can be assessed as the computation is carried out, or, some approximate 
formulae must be developed which are valid within as wide a range of the para- 
meters as possible. In the application to a specific example, the latter method 
would be used, while the former method could be resorted to as a check in cases of 
doubt, if alternative exact theories are not available. 


One point of interest may be noted. It is that for wings with subsonic leading 
edges the accuracy of any results at Mach numbers very near to unity, obtained by 
an Evvard technique, is questionable. However, since it is known that the exact 
linearised solution at M=1 is given by the slender-body theory, it should be possible 
to interpolate between M=1 and the Mach number range in which Evvard’s 
method can be applied. It should be emphasised that the validity of any results 
obtained by this method are bounded by the conditions and assumptions on which 
the linearised equations are based. Although the linearised equations do not 
represent the physical flow conditions at M=1, the procedure affords a necessary 
mathematical check of results in this transonic range. 


NOTATION 


Cartesian co-ordinates 


u,v characteristic co-ordinates 
p(x, y) pressure distribution 
Mach number 
speed of sound 


a 
(“) =f ‘e symmetric equations of wing leading edge 


SwnsSpn areas of integration 

a(u,v) incidence 

A(u,v) upwash velocities 

@(u,v) velocity potential 
B = (M?-1)! 


6 angle of sweep. 


2. Mathematical Formulation 


On the application of the linearised theory to problems of supersonic flow 
around finite thin wings, it is known that the value of the perturbation velocity 
potential ¢ (x, y, z) on the upper surface of the lifting discontinuity, taken on the 
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plane z= +0 (defining the supersonic flow of velocity V to be in the positive 
x-direction), is given by 


A(E, 0) d&dy 
(x, y,+0)= | (1) 


Ss 


where 


(i) | S is the area enclosed by the intersection of the fore Mach cone from the 
point (x, y, 0) with the xy-plane; 


(ii) A(E,) is the downwash velocity on the xy-plane and, in particular, on 
the wing surface itself A(¢,») is given directly by the local angle of 
incidence. 


On transformation of the Cartesian co-ordinates (€,7) to characteristic 
co-ordinates (u, v), defined by 


M ,. 
u= (€ — Bn), v= 2p Pu) 


B 1 @) 
or M (u+v), & u), 
equation (1) reduces to 


with the transformed reference point (u, v) corresponding to (x, y, +0). 


Fig. 2. 
Notation. V 


v=f(u) 
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This paper is concerned with thin aerofoils symmetric about the x-axis (i.e. 
zero yaw) with subsonic leading edges and pointed noses. Hence, the area of 
integration S$ is given by OX u=<=u,0<v<vz, as shown in Fig. 2, and the 
equations of the leading edges are denoted by 


Writing S=Sy+Sp,+Sn2, 


as shown in Fig. 2, where S, is the contribution to S of the wing surface, Sp, and 
Sp. are the areas in the neighbourhood of the leading edges, equation (3) becomes 


SpitSpo 


where 2 (u, v) is the local angle of incidence on the wing surface. 


It is known, however, that this integral can be reduced further by Evvard’s 
technique to the form 


V2(u, v) dudv V2(u, v) dudv 
| “w2 
A, (u, v) dudv A, (u,v) 


D1 “D2 


where the areas S,,. Sw., Sp, and Sp. are shown in Fig. 3; and A, (u,v), A. (u,v) 
are the upwash velocities in the areas S),, Sp. respectively. 


V Fig. 3. 
Areas of integration. 
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The development of equation (5) can be followed in two courses, as discussed 
in Section 1. First, one can attempt to carry Evvard’s method to the next stage, 
reducing the integrals over the upwash regions to integrals over the wing surface, 
which in general will have to be integrated numerically. Alternatively, estimations 
can be made of the upwash velocities A, (u, v) and A, (u, v) for the direct integration 
of equation (5). 


3. Estimation of Upwash Velocities 


To determine the upwash velocities A, (u,v) and A, (u,v), equation (1) is 
applied to a point off the wing surface but within the bow Mach cone. This 
configuration is indicated in Fig. 4. 


By reference to Fig. 4, equation (1) becomes 


Sw pit Spe 


For points (up, vp) off the wing surface, 
(Up, Vp) =0 ‘ ‘ (7) 


since the velocity potential is continuous in this region and anti-symmetrical with 
respect to the uv-plane. 


vif (u) u= f (v) 
vef'w 


Fig. 4. 
Notation for a point off the wing surface. 
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Therefore, combining equations (6) and (7), 


= v=f(u) D 
(u, A,(u, v) dv +f A, (u, v) (8) 
(vp - 7) (up (vp - 


in the notation previously introduced, using the idea of an operator f denoting the 
leading edge equation, so that u=f (v) is equivalent to v=f~' (u). 


Thus, since up is a general point on the upwash field Sp., the integrand with 
respect to u in equation (10) may be put equal to zero; hence 


v= v=f(u) 


(u,v) dv (u,v) dv (u,v) dv (9) 
v=/(u) v=0 liu) 


A similar expression can be written for the area Sp,, giving, with equation (9), 
a pair of simultaneous equations for A, (u,v) and A, (u,v); but this paper is only 
concerned with a wing which is symmetrical about the line «=v, and the relation- 
ship between the two upwash velocities is 


If equation (9) is rewritten in the form 


(u) v=f(u) 
A, (u,v) dv V2 (u,v) dv A, (u,v) di 


v=f(u) v=0 


then this could be regarded as an Abel type of integral equation for A, (u,v). It is 
known that if 


x 


Hence. the solution of equation (11) is 


=i (u) 
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This equation can be regarded as an integral equation for A, (u, v) by using the 
symmetrical relationship of equation (10). The second term depending on X, (u,v), 
however, is in general smaller than the first term, depending on the wing incidence, 
for Mach numbers greater than unity, and may be neglected for a first approxima- 
tion to the upwash problem. Further approximations may be obtained by a 
successive approximation technique of substitution of the first value of A, (u,v) in 
the second integral. 


4. Method 1 


In this method we attempt to extend Evvard’s idea of replacing the upwash 
integrals by integrals over the wing surface. We begin by substituting equation (12) 
for A.(u,v), and the similar symmetric integral expression for A, (u,v), in 
equation (5). Consider first the integral 


D2 Niu) 


where u,=f (v,)=f?(u), shown in Fig. 3. Then, on substitution of equation (12) 
in equation (13), the expression (14) is obtained, on integration with respect to v, 
and then changing the variable v’ to v. 


us f(v) 


Fig. 5. 


Areas of integration. 
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=f(u) 


v=f(u) 


Therefore, using equation (14), equation (5) reduces to 


Va (u,v) dudv V2 (u,v) dudv 
(u, ) <M = [(u—u) (v7 | | {(u—u) (v —v)]? 


“wi 


w3 


ws 


Spy 
(15) 


where the areas of integration are shown in Fig. 5. 


Theoretically it is possible to repeat the process of substitution for A, (u, v) and 
A, (u, v) in equation (15) and so reduce further the integrals of the upwash velocities 
to regions of the wing surface. A reasonable approximation would result by 
neglecting the integrals over Sp, and S»., giving a value slightly greater than the 
accurate solution, but the integrations over the remaining areas would in general 
have to be computed, preferably in the manner suggested by Etkin and 
Woodward™. A further approximation could be developed by replacing A, (u, v) 
and A, (u, v) by the first term in equation (12). 


As pointed out in the Introduction, Evvard’s assumption involved the integral 
over S,,,, while Etkin and Woodward included also the integral over S,.. However, 
for more reliable results recourse must be made to an equation similar to (15). 
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Fig. 6. 
Areas of integration. 


The question of accuracy arises if the process indicated by equation (15) is 
stopped at any particular stage and, in general, the maximum error will appear on 
approaching M=1. As pointed out in the Introduction (Section 1) a mathematical 
check is available at M=1 from slender-body theory, so that all that is necessary is 
to ensure that the integrations give solutions as true as possible for conditions as 
near to M=1 as possible. In general, one could only do this by actually computing 
specific examples, but it is interesting to discuss the convergence of a similar 
expansion to equation (15). It is assumed that the velocity potential is given by 


1 Va(u, v) dudv 
9. —M ) 
-SwotS3- Sy 


where Sy,, Swe, S;, etc., are indicated in Fig. 6. 


For a triangular flat plate, z=constant and v=ku, u=kv are the equations of 
the leading edges (k being a constant depending on the angle of sweepback and the 
Mach number); hence 


The pressure loading then becomes 


(a, 2pv =2 pv (#4) + 2) 
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At M=1, since k —> 1 as 8B —> 0, so that (1 - k)/8 is finite in the limit, the pressure 
loading on the centre chord line u=7 is 


p (u,u)=K {2-0:34+0-:20— .... +(—1)"2 } 


where K depends on pa’, and the angle of sweep. Although the convergence of 
this series is slow the first three terms give an answer to within four per cent. 


On the basis of this example it could be argued that the approximation to 
equation (15), as suggested, forms the starting point for calculating pressure loadings 
for all practical purposes. 


5. Method 2 


In Section 4 a rigorous approach to the problem is presented which involves 
rather complicated integrals, leading to a numerical approach which can be used to 
give any desired accuracy. However, in general, this would be too cumbersome, 
so that an alternative, if more approximate, technique should be developed. This 
possibility is investigated in this section. 


The method is to obtain expressions for the upwash velocities directly from 
equation (12), substitute them in equation (5) and integrate approximately. 


Now equation (14) states that 


+an integral depending on A, (u. v). 
It is possible to neglect the effect of the second integral depending on A, (u, v), 
but the percentage effect can be calculated by a successive approximation technique, 


as suggested previously. The integral over the wing can be worked out fairly easily; 
for example 


(i) flat plate, i.e. z=constant 


(ii) linear camber, i.e. z=C,u+C.v, then 


Xr, (u, y= 


v-f-* (u) f-* (u) 


_ 2VC, WP 
3x v—f-*(u) 


and similarly for more complex cambers. 
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The second approximation is found by substituting the first approximation 
derived already in the second integral in equation (12); hence for a flat plate wing 


For a flat plate with straight edges 
(")=1(:) 
(49. 


where k is a constant and equal to (1—f tan 6)/(1+3 tan 6), @ being the semi-apex 
angle. Equation (18) reduces in this case to 


A. (u, v)= { ~tan-! E (i 


kv -u kv -u 
u (1 -(1—k’*)*] }! 


The maximum effect of the third term compared with the first term is for u= kv, 
and then the ratio of these terms is 


For k=073 0-025, 
k=0°5, = 0-085, 
k=0-7 0-30, 


so that it seems that, for k=0-5, the third term could be neglected in comparison 
with the first term, but for higher values of k it should be taken into consideration. 
It is worth noting that k=0-5 corresponds to @ tan @=4, so for a 45° sweep this 
corresponds to Mach 1-05; hence the range of k greater than 0-5 only contributes to 
the range of Mach number 1-00 to 1-05 for this particular wing. An approximate 
theory which is valid down to a Mach number near to unity is sufficient for most 
practical cases. For the rest of the analysis this third correction term in 
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equation (19) will be neglected. The foregoing example is taken as a test case in 
respect that other problems, such as cambered wings, curved leading edges, and so 
on result in terms of the same relative order in an equation corresponding to (19). 


It follows from the foregoing argument that in equation (15) the terms involving 
the upwash velocities can be neglected to the same degree of approximation as 
intimated. 


6. Determination of Pressure Loading 


The integral for the velocity potential is given in equation (5), and may be 
written 


1 
_a(u,v) dudv 
where h=V | | (21) 
Swi-Swo 
A (u, v) dudv 
Spi 


The areas Sy,, Sw, Sp, and Sp, are shown in Fig. 3. The wing incidence =z (u, v) is 
known for any given wing. The upwash velocity is determined by the method 
suggested in Section 5. 


The pressure loading p (u, v) is given by 
(u, v)=2p V 


on account of transformation equation (2). This can be applied to equation (20) 
in two ways. Either the integration is performed first, as defined by equation (20), 
followed by the operational equation (23) for each specific example, or, general 
expressions can be developed for p (u,v) by differentiating equation (20) under the 
integral signs. In both cases the expression involving /, (equation (21)) will be 
composed of standard integrals for most leading edge contours and cambers. The 
evaluation of equation (23) is complex but a slight modification may be introduced 
which reduces the integrals to a more amenable form. Consider, for example, the 
integration over area Sp, (see Fig. 3). Since the variation of u is within 
the range OSu<f(v,), or OSu<f?(¥), and f?(u)<<u, the function 
u—u may be replaced by u-f?(u) and taken outside the integral sign. This 
approximation may be regarded as being of the same order of error as those 
incurred in the calculation of the upwash velocities A (u, v). 


These ideas have been applied to the calculation of pressure loadings over 
particular wings. 
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6. Triangular Flat Plate 
6.1. THE GENERAL CASE 


In this case the equations of the leading edges are given by 


_ 1-Atané 


Btand’ being the semi-apex angle 


where, as previously noted, k 
Mach | corresponds to k= 1, and occurs when the leading edge becomes sonic. The 
formulae of the pressure loading are simple in this case. 


The graph of C,/(22 tan 4) against 8 tan 6, where C, is the pressure coefficient 
on the centre chord line, is shown in Fig. 7, derived by the foregoing method. It is 
compared with the exact solution, and also with that of Etkin and Woodward. 


6.2. CRopPED DELTA FLAT PLATE 


The application of this approximate theory to the problem of the severely 
cropped delta wing has greater practical significance than that for the triangular 
wing, since the exact solution of the problem involves the superposition of conical 
fields. The wing contour is defined by 


“ =k, for the leading edges 


_ (u\ v+k, 
and = (“) «+, for the wing tips, 


together with the Mach regions of influence A, B and C, as shown. 


veu+k, 
/ 
8 
\ 
\ / 
\ / 
V 


The pressure distribution in area A is given by the solution for the triangular 
wing in Section 6.1. The pressure at a point (u,v) in areas B and C follow directly 
from this solution, since the related point (u,,v,) lies in area A, so that the upwash 
velocities around the leading edge only are needed. 


The loading has been calculated for a wing of 45° sweep, unit tip chord, root 
chord=2 at M=1-1, and is shown in Fig. 8. 
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6.3. FLAT PLATE WITH CURVED LEADING EDGES 


In this case the equation of the leading edge will be of the form 
v=f(u)=k, 


where the k,’s are constants depending on the wing contour and Mach number. It 
can be shown that the integrals reduce to standard types, especially if the wing 
contour changes slowly so that v =f (u) can be expressed by the first few terms of its 
Taylor series with reasonable accuracy. 


6.4. CAMBERED AEROFOILS 


The procedure for cambered aerofoils is exactly the same as that for the flat 
plate aerofoils. The approximations intimated are sufficient to reduce the integrals 
to standard forms if the incidence distribution is given as a series, although the 
number of integrals to be evaluated rises appreciably as the complexity of the 
incidence increases. However, the method is straightforward, even if tedious, and 
has the fundamental advantage of being completely general. 
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Non-Linear Elasticity in Aluminium Alloys 
A. B. WATTS*, B.Sc.(Eng.), Ph.D. 


SUMMARY: Tests have been carried out on round bar specimens of various 
strong aluminium alloys (unclad) in the fully heat-treated condition in both 
tensile and compressive loading. A careful study of the results establishes 
beyond reasonable doubt that there is a progressive change in the tangent 
elastic modulus with stress. The variation is appreciably greater than that 
predicted by the consideration of second-order effects in the conventional stress 
and strain representation. It is found that the modulus decreases with tension 
and increases with compression by about 5 per cent at 30 tons/in.?, the form 
of the variation being a continuous function. Tests on two alloys of aluminium 
other than those of the normal zinc-bearing and the copper-bearing types were 
also carried out. The first of these was a simple fully heat-treated binary 
alloy of 4 per cent copper which was very similar in behaviour to D.T.D.364; 
the second was a 5 per cent alloy (N.S.6-4+H) which, after its initial loading, 
showed only very slight elastic modulus variation in the subsequent loading 
and unloading lines. An additional test to compare naturally and artificially 
aged D.T.D.364 did not reveal any difference in the degree of modulus variation. 


In view of the present experimental findings, the use of the “second” 
modulus line in the proof stress determination of unclad strong alloys should 
be adopted generally, if practicable. A high level of internal stress can cause 
slight errors in this method because the elastic limit is then prematurely reached. 
As an occasional check on the method when applied to routine mechanical 
testing, the elastic component of strain can be found for the last point on the 
loading line by returning to the zeroing load. 


In the Appendix the difference between true and nominal stress-strain 
relationships is calculated. The meaning and theory of modulus variation is 
also discussed and aluminium alloys are compared with a high-tensile steel. 


A simple expression for the form of the variation is proposed. 


1. Introduction 


The non-linear elasticity of metallic materials is a direct contradiction of the 
simple laws which are commonly regarded to hold true. In attempting to decide 
for any material whether or not there is non-linear elastic behaviour, it is important 
to eliminate a number of factors which can, and do, modify the shape of a stress- 
strain curve such as is obtained in a proof stress determination. If a particular 
material does exhibit a genuine curvature of the purely elastic portion of its stress- 
strain curve when all these factors have been eliminated, then the result has 
far-reaching implications in several different fields. 


*The work described in this paper was done while the author was a research engineer with The 
British Aluminium Company. 
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6.3. FLAT PLATE WITH CURVED LEADING EDGES 
In this case the equation of the leading edge will be of the form 
where the k,,’s are constants depending on the wing contour and Mach number. It 
can be shown that the integrals reduce to standard types, especially if the wing 


contour changes slowly so that v =f (u) can be expressed by the first few terms of its 
Taylor series with reasonable accuracy. 


6.4. CAMBERED AEROFOILS 


The procedure for cambered aerofoils is exactly the same as that for the flat 
plate aerofoils. The approximations intimated are sufficient to reduce the integrals 
to standard forms if the incidence distribution is given as a series, although the 
number of integrals to be evaluated rises appreciably as the complexity of the 
incidence increases. However, the method is straightforward, even if tedious, and 
has the fundamental advantage of being completely general. 
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